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SUMMARY 


A complete second-order closure program for the investiga- 
tion of compressible, turbulent,^ reacting shear layers has been 

developed.' \T.he; e qua tio ns.; ' 'for _ the ..means. . .and . i K&^sejcond-__ .• XT’ u 

order correlations are derived from the time-averaged Navier-iJ 
Stokes equations. The equations contain third-order and higher- 
order correlations, which have to be modeled in terms of the 
lower-order correlations to close the system of equations. The 
fluid mechanical turbulence models and parameters are the same 
ones used in previous studies of a variety of incompressible and 
compressible shear flows. A number of additional scalar correla- 
tions have to be modeled for chemically reacting flows, and a 
"typical eddy" model has been developed for the joint probability 
density function for all the scalars. A simpler model in which 
the higher-order scalar correlations are set zero is also being 
used/ The program is capable of handling multi-species, multi- 
step chemical reactions. 

Computations .have been carried out using the second-order 
closure reacting shear layer (RSL) program to study nonreacting 
and reacting flows. The calculations for a hydrogen-air diffu- 
sion flame demonstrate the importance of the inclusion of the 
concentration fluctuation correlations in the chemical source 
terms. The neglect of the correlations leads to the prediction 
of a very thin' reaction zone^ whereas, when the unmixedness 
effects are included in the analysis, the correct turbulent flame 
structure is calculated. 

The concept of the "typical eddy" model has been verified 
for low heat release reacting flows, and substantial progress 
has been made towards the construction of a complete model in- 
cluding the density correlations. 


INTRODUCTION 


Predicting the characteristics of turbulent flames (both 
diffusion flames and premixed flames) is a problem of major im- 
portance in many fields of engineering and has been studied off 


and on for almost half a century. The problem basically remains 
unsolved, though a'-considerable amount of engineering progress 
has been achieved by experimentation and empiricism. However, 
a truly predictive solution of the turbulent reacting flow 
problems has not yet been achieved. Actually, even the simpler 
problem of predicting turbulent nonreacting flows has not been 
solved, though considerable advances have been made in recent 
years in the development of higher-order closure models. The 
use of these new closure models for nonreacting flows have 
demonstrated improved predictive capability compared to the eddy 
viscosity or first-order closure models. A.R.A.P. has been one 
of the active groups in the development of multi-equation, com- 
plete second-order closure modeling of turbulent flows. The 
application of these higher-order closure procedures to turbulent 
reacting flows is a reasonable extension of current technology. 

This report describes the recent advances made at A.R.A.P. on 
the methodology for the study of chemical reactions in turbulent 
media. 

Turbulent flows involving chemical reactions occur in many 
situations including industrial and home furnaces, chemical 
process plants, various propulsive devices such as jet engines, 
rocket motor-ram jets, chemical and gas dynamic lasers, and 
exhaust plumes and wakes of reentry vehicles. In effect, al- 
most every practical flow system involving chemical reactions 
involves turbulent flowfields. The interaction 'between the 
chemistry and the turbulence is of significant importance in 
many of these systems. The interaction is of importance in 
determining combustion efficiency and pollutant formation as 
well as many other combustion characteristics, such as igni- 
tion and extinction behavior, flammability limits, combustion 
stability, combustion noise, etc. The dual and often contradictory 
objectives of increasing the combustion efficiency and controlling 
pollutant formation require a more fundamental understanding of 
the Combustion process; and one of the major features . that has : 
been ignored in most previous studies is the effect of the 
turbulent fluctuations of concentrations and temperatures on 
the chemical reactions. 


The main reason this problem has been ignored in the past 
is due to its complexity. The need for understanding the coupling 
between turbulence and the chemical processes has been appreciated 
and acknowledged by the early combustion researchers. However, 
until very recently the analytical and computational capabilities 
as well as the diagonistic tools to adequately characterize 
turbulence in reacting flows were not available. Considerable 
advances have been made in recent years in these areas and it 
now appears feasible to:' develop a rational model for turbulent 
reacting flows using the second-order closure approach. 
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In a turbulent diffusion flame, the reaction between 
initially unmixed reactants (i.e., species a and g ) takes 
place primarily at the edges of the streams of a and g that 
are twisted together and intertwined as a result of the turbulent 
mixing process. In this situation, although the time average 
concentrations of the species a and g at a point in the flow 
may be quite substantial, the reactants may still not be mixed 
on the molecular level. The local rate of the chemical reaction 
is governed by the degree of molecular mixing, and is usually 
significantly lower than the rate that would be calculated using 
the time average concentrations of the species at that point. 

The effect of this molecular unmixedness on the reaction can be 
incorp orate d by considering the effect of the species correlations 
like a ’ g * , on the mean chemical source terms. Terms of. this 
kind are calculated in a second-order closure analysis. 

A number of efforts are now underway to develop 'a' * '• ' .. 

model of this kind. Donal dson and Hilst ( r £ f . ' 1A studi ed, 1- 

low heat release chemical reactions and demonstrated the im- 
portance of keeping track of the species correlations. An 
explicit closure model for third-order species correlations was 
developed, but this is only valid for isothermal flows. Spalding 
(ref. 2) formulated an eddy break-up model to study the effect 
of turbulence on the overall combustion process, but the model 
is only valid for large turbulence Reynolds numbers and ignores 
the chemical kinetics effects. The major new problem in the 
analysis of turbulent reacting flows is the modeling of the higher 
order scalar correlations that arise from the chemical kinetic 
source terms. There seems to be a general acceptance that the 
desirable procedure is to model or calculate the probability 
density function (pdf) for the scalars. This procedure has been 
adopted by A.R.A.P., as well as by Rhodes et al. (ref. 3), Bray 
and Moss (ref. 4), Libby (ref. 5), Lockwood et al. (ref. 6), 

Borghi (ref. 7), and Spalding (ref. 8). 

Most of these approaches still need to make many simplifying 
assumptions to make the problem tractable. Bray and Libby use 
the fast chemistry assumption and deal only with a one-dimensional 
pdf. The approach of Borghi is very similar to ours for the case 
of a one-step kinetic reaction. However, he also replaces the 
joint pdf of species and temperature by a one-dimensional pdf 
which is a very restrictive assumption. The one-dimensional pdf 
being popularly used in most of these studies is a clipped 
Gaussian with delta functions at the two edges. Spalding has 
recently proposed a Lagrangian theory of turbulent combustion 
which incorporates the kinetic effects of multi-step chemistry 
and involves calculating the pdf at various points by following 
the history of various eddies. 

The A.R.A.P. approach models the joint pdf of all the scalars 
in turbulent reacting flows. A procedure has been developed for 
constructing the pdf using the available information in a second- 
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order closure analysis. The model is called the "typical eddy" 
model and involves representing the pdf by a set of delta func- 
tions of variable strengths and positions in the scalar phase 
space. The strengths and positions of these delta functions at 
every point in the flow are determined from the predicted values 
of the means and second-order correlations at that point. 

Kewley (ref. 9 ) has recently used our concept of the "typical 
eddy" model in a study of turbulence effects in a chemical laser. 
The RSL program has also been used for an extensive study of tur- 
bulence effects on DP chemical laser flowfields which is detailed 
in Reference 10. 


This report describes the models and the A.R.A.P. Reacting 
Shear Layer (RSL) program code. Results of model and program 
verification studies and calculations for the hydrogen-air dif- 
fusion flame at conditions corresponding to the experiments of 
Kent and Bilger (ref. 11) are presented. The calculations 
demonstrate the importance of including the scalar correlations 
of the species and density in the chemical source terms for 
finite rate chemistry calculations to avoid significant errors 
in predicting the structure of turbulent flames. 

SYMBOLS 


a 

A 

b 

D 

g 

h 

k 


Pj j ^2 > ^3 


r l/2 

r* 

Re 


turbulence model parameter 

turbulence model parameter 

turbulence model parameter 

diameter 

matrix tensor 

sensible enthalpy 

thermal conductivity 
reaction rate 

pressure 

turbulence model parameters 
turbulence kinetic energy 
radial coordinate 

half radius of mean velocity profile 

2 12 

radius for which q = tt 9 

n A unax 

Reynolds number 
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s , Sip 
T 

u i or u, v, w 

V 

c 

x,X 

y 

a,8,Y 

'ei 

1 

0 

A 

X 

y y y * 

V 

p 

Superscripts 
", < > 


turbulence model parameter 
temperature 
velocity components 
turbulence model parameter 
axial coordinate 
normal coordinate 

f 

species mass fractions 

cell size in "typical eddy" model 

T - T 

ambient 

turbulent macroscale 
turbulent microscale 
molecular viscosity coefficients 
kinematic viscosity = y/p 
density 

denotes time average 

denotes fluctuation about the mean value 


SECOND-ORDER CLOSURE MODEL AND PROGRAM DEVELOPMENT 


Basic Equations 


The general conservation equations for a turbulent, com- 
pressible, multi-component reacting flow system are presented in 
Appendix A. These are the Navier-Stokes equations for overall 
mass continuity, conservation of momentum and sensible enthalpy 
and for species continuity. The perfect gas equation of state is 
used. The following assumptions are used in the derivation of 
these equations: 
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1. The gas is a continuum. 

2. No body forces are present. 

3. Radiant heat transfer is neglected. 

4. There are no overall mass sources. 


5. Mass diffusion is only due to a concentration gradient 
Pick's diffusion law is used. Thermal and pressure 
diffusion are neglected. 

6. Each component of the gas mixture is thermally perfect 

The dependent variables in these equations are decomposed 
into a mean and a fluctuating part, and equations for the mean 


quantities u. 


a 


0 


and various second-order 


uTTP 


u! a ' 


oP r T r 


etc . are 


correlations, such as uTuT , , 

1 i j i i 

derived. This system of equations is not complete as the equa- 
tions contain a number of third-order and higher-order correla- 
tions that have to be modeled to obtain a closed set of equations. 
This is the familiar closure problem in the analysis of turbulent 
flows. A second-order closure approach involves modeling these 
higher-order correlations in terms of the means and the second- 
order correlations. The basic task in accomplishing a second- 
order closure solution of turbulent flows is the development 
of suitable models that remain invariant with respect to changes 
in the flow geometry. The general model parameters can then be 
obtained by experiments on simple, standard flowfields and the 
model can be used as a predictive tool for new and more complex 
flow problems. 


The following additional simplifying assumptions are used 
in the present study: 

1. The mean flow is steady. 

2. The mean pressure is constant in the direction normal 
to the mean flow. 


3. Fluctuations in u 


k are related to the 


temperature fluctuations by the expressions of the form 


' V T T ' + l y a a ’ - 

a v •- 


+ y in. 

1 3a 


a 


a 


a 


V ',i = U T T |i + I “a “li 
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k are 


4. Third-order correlations involving y , y* 
neglected. 

.. The various models!' for. 'higher-order correlations that 
have been developed are reviewed- .;in "‘the next' --section . 


Second-Order Closure Models - Nonreacting Plows 


Considerable advances have been made in recent years in the 
development of second-order closure models for nonreacting 
turbulent flows, by many groups of investigators (refs. 12,13,14). 
Donaldson and his colleagues at A.R.A.P. have been involved in 
the development of a multi-equation complete second-order closure 
procedure for turbulent incompressible and compressible flows 
(refs. 15,16,17). The models are selected in the simplest 
possible way consistent with tensor invariance and available 
experimental data. The model must remain invariant with respect 
to changes in the flow geometry. The invariant modeling is an 
important feature in the sense that our goal is a model which, 
although it is semi-empirical, has no varying constants which 
must be determined for each new flow. The model must be 
applicable to a wide class of fluid flows. 

The models for various classes of terms that appear in the 
second-order closure equations for nonreacting flows are briefly 
described in this section. The same identical models and model 
constants are also used for these terms in problems of reacting 
flowfields. The models and the philosophy behind them have 
v been previously described in greater detail in a number of reports 
and publications (refs. 17,18,19) . . For the sake of completeness 
of this report, some of this information has also been included 
here . 


Modeling of pressure fluctuation correlations .- The correla- 
tions involving pressure fluctuations are among the more difficult 
terms to model in a second-order closure analysis. The results 
for compressible free shear layer flows presented in Reference 15 
indicated the need for improved modeling of the pressure diffu- 
sion correlations. New models for various pressure correlations 
in compressible flows have been recently developed. The models 
are obtained by following a consistent procedure beginning with 
the Poisson equation for the pressure, 


a 2„ 

i_E. = 

8x? 




9u i\ 

8 x j/ 
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The equation is resolved, into its mean and fluctuating parts 
to obtain an equation for p' . This can be formally integrated 
to obtain p' and then expressions for various correlations can 
be written. The integral cannot be solved, but with the use of 
certain approximations, it can be used to suggest the form of 
the modeled terms. The procedure is described in Reference 20. 
The models are incorporated in the program code but currently 
are not being used in the computations in the form given below, 
as parts of the model have not yet been . completely tested. Many 
combustion problems involve relatively low speed flows, and the 
program is applicable to the study of these systems. However, 
the extension of the program to high Mach number flows requires 
the validation of the constants to be used with portions of the 
pressure models given below. 


The pressure fluctuation correlations arise from the pres- 
sure gradient terms in the momentum and energy equa tion s (Eqs. 
A-3 and A-4)*. The term is typically of the form c|)'p' k 


*'p; k - (p'*' ) > k - P'4»; k 


Models have been developed 
hand side. The terms appear in 

equations. For example, p'u! 

J 

tions for u ! u ' (Eq. A-16) and 
J 

p'h' and p'h' . appear in the 
(Eq. A-17). 


for the two terms on the right- 
many of the second-order closure 

and p'u! k appear in the equa- 

upT*" (Eq. A-18). Similarly 

equations for u'.h' and h'h' 

J 


Pressure diffusion terms 


p'u| = P-j^pA 


- « 2 m 


3x. 


n L 


i(u’V) 


\ { 6 m (u i u,P) , P + g mi ( 


U 'V%} 


+ I { 6 m u i u 'p - «W “'"“Ip } + p 2 PAq(uf“’ P ). >p 


(1) 


p'h' = P, pA 


~ 3u 

-.2 m 

9x n 


(h'u' n ) - 6 n (h'u' 5 ') 


,m 3 


+ i 6 n h'u' 


3 m 


+ P 2 pAq(h'u' ) % 


( 2 ) 


*The A in the Equation numbers refers to Appendix A. 
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Tendency-towards-isotropy terms 


p ’ u i,j 


_ - m 

, n . 

1 

,n 

P 3 P 3x 
J n 

®mj u U i - 

3 

6 

m 

1-- ~ 

n , . 1 r n 


2 

3 6 U U 

u ; + 9 6 m 

S ij 

q 


u : u : 
i j 


1 23. 

2 A 


i J 


3* s ij q 


(3) 


3u 


p ' h ;i = p 3 p 


m 


n 


T u ' h 
3 m i 


j 6? u'h' 
2 1 m 


I s mi u 




- 0.8 £9- ujh' 


(4) 


In these models, A is the turbulent macro-scale length. 
The major change compared to the models used in References 15 
and 18 is the addition of terms involving the mean strain 

3u /3x . The constant P~ has been selected to be 0.1 from 
m n 2 

the studies in Reference 18. The constants P-^ and P^ 

associated with. the mean strain terms are still under investiga- 
tion and have to be calibrated using data for compressible shear 
flows . 


Modeling of third-order velocity correlations . -These terms 

represent the process by which various second-order correlations 

are transferred from one part of the flow to another without any 

net production or loss. These terms appear in all the secon d- 

order correlation transport equations. For example, ulu'.u.' 

i 3 ^ 

and u|u'. p' appear in Eq. A-l6, u ! h ' h' and u!p'h' are present 

in Eq.A-17, etc. The most popular "‘'modeling of £hese transport 

terms is as a gradient diffusion process. A number of different 

gradient diffusion forms has been proposed (ref. 17). Following 

previous studies (refs. 15,18), a simple gradient diffusion model 

is adopted that satisfies the tensor symmetry, and has a scalar 

diffusion coefficient. The models are of the form: 


u ! u !u/ 

l J k 


U i U j*' = 


u^<|> ' 4> ' 


- v c Aq [ (u i u j ) .k 

-v o Aq [aqrn^ + 
(FT 1- ) J 


+ (u k u i\j 


-V Aq 
c 


(5) 

( 6 ) 
(7) 
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In Eqs . (6) and (7), <j> * represents any scalar variable 

(T' , h r , c^ , p' , etc.). V c was established equal to 0.1 by 
comparing model results with experimental data in free jets and 
shear layers. Due to the lack of detailed experimental data in 
compressible and reacting flows, the same diffusion constant is used 
as an initial estimate for the velocity diffusion of all correlations 


An alternate set of model constants for pressure and velocity 
diffusion has been used in atmospheric flow studies. In these 
studies (refs. 16,17), P 2 = -0.3 and V c = +0.3 • There are now 
some indications that this is a somewhat better model on the 
basis of some recent studies of mixing layer flows of two species, 
but further ‘work is still continuing to settle this difference. 

Dissipation correlations .- The dissipation correlations of 

the form g mn (uT! u7 ) appear in the equation for u u,' (Eq. 

0 j , m k,n ^ p j k ^ 

A-l6) when the viscous terms u,' x 0 and u [ 0 are expanded. 

These terms measure the effect of viscous decay on the second- 
order correlations. Viscous dissipation is expected to be the 
major loss mechanism for the turbulence kinetic energy and other 
correlations. Lewellen et al. (ref. 21) suggest that an isotropic 
dissipation model is more suitable than the completely anisotropic 
dissipation model used by Donaldson (ref. 18). However, it is 
likely that low Reynolds number dissipation will be anisotropic 
and the following form of the model, which goes to the correct 
limits for small and large Reynolds numbers is now used. 


g mn (in 


UJ 


j k,n 


) 


g 


Jk 


3A £ 


A' 


( uj u k 


- g 


jk 


2 

a_ 

3 


( 8 ) 


X is the Taylor microscale length, and Rotta's form for it is 
used. 


A 2 = 


a + 


bqA 

v 


Detailed studies of the flat plate 
layer (ref. 22) have led to a revision 
2.5 to 3.25 . The constants a and b 


incompressible boundary 
in the value of a form 
are : 


a = 3.25 
b = 0.125 

Dissipation correlations of the same type appear in all the 
second-order correlations from the molecular transport terms. 

The models for all the other dissipation correlations involving 
scalar fluctuations are of the following form: 


mn rp ! m t 
,n 


rp ! rp t 



(9) 


s = 1 . 8 
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( 10 ) 


mn 

g 


! T 
i,m , 


u! T ' 



s T = 1.0 


The models for the dissipation of the other scalars are 
of the same form as Eqs. (9) and (10). 


Turbulent macro-scale A .- It is necessary to specify the 
turbulent scale A , used in many of the modeled terms, in 
order to complete the closure of the second-order correlation 
equations. For nonreacting flows, both an algebraic specifica- 
tion of the scale based on the gross features of the flow- 
field, and the use of a modeled dynamic differential equation 
to calculate the local scale at each point in the flowfield have 
been used. The latter procedure does lead to somewhat better 
results, but up to this time it has not been possible to develop 
a universal model scale equation. Some of the parameters in the 
scale equation have to be changed for different flow geometries. 
This problem can be readily appreciated when one realizes that 
all of the terms in the scale equation have to be modeled. The 
scale equation contains much more arbitrariness than the other 
second-order correlation equations. 

With the uncertainty involving the use of an incompressible 
dynamic scale equation and the lack of a scale equation formula- 
tion for compressible and/or reacting flows, it was decided to 
simply use an algebraic specification of the macroscale in the 
reacting shear layer (RSL) computer code. The scale is determined 
using either the mean axial velocity, u , profile or the profile 

2 

of the turbulence kinetic energy, q , For axisymmetric flows, 

1. A = 0.5 r ]L /2 3 where r i/2 is radius of 

the mean axial velocity u . 

2 12 

2. A = 0.2 r* , where r* is the radius where q = K q . 

*4 ITlcLX 

The use of the turbulence kinetic energy profile for 
determining A gives better results for nonreacting flows. At 
the present time there is no clear conclusion as to the preferable 
scale procedure for reacting flows, due to uncertainties in 
other modeling approximations . The scale based on the mean axial 
velocity profile has been used in most of the current reacting 
flow calculations. 


Second-Order Closure Models - Scalar Correlations 

in Reacting Flows 


The major new difficulty in the second-order closure com- 
putation of turbulent reacting flows over that for nonreacting 
flows is the need for modeling of scalar correlations, such as 

k’a' , k ' a ' 3 ' , a 7 a ' 3 ' , etc. that appear in the equations. 
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The "typical eddy" model for the joint probability density 
function of all the scalars ( p , h , T , ou , k ) has been 
developed by Donaldson and Varma (refs. 23,24) that appears 
to be very promising for the modeling of all such scalar 
correlations. 

Consider a simple one-step reaction involving two reactant 
species a and 3 and one product species y , i.e., 
a + 3 y . Even for this very simple chemical system, a very 
large number of higher-order scalar correlations appear in the 
equations and particularly in the chemical source term. Instead 
of attempting to individually model all these terms, there is 
a general acceptance among many of the research groups active 
in this area, to model or calculate directly the probability 
density function (pdf) of the scalars. In general, the pdf will 
be a continuous function in the scalar phase space and attempts 
have been made to model the pdf by various continuous functions 
(refs. 4, 6, 7) for the case of one-dimensional pdf's. However, 
for a multi-dimensional joint scalar pdf it is much more diffi- 
cult to guess the appropriate shapes for the pdf. The "typical 
eddy" model proposes to represent the joint pdf by a set of 
delta functions of variable strengths and positions in the 
scalar phase space. The sketch illustrates how the important 
features of a complicated continuous pdf may be represented by 
a series of delta functions, and also shows an alternate "box" 
representation for the two-dimensional delta function pdf. 

Once the joint pdf has been modeled, it is possible to calculate 
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all the required higher-order scalar correlations and to close 
the set of second-order correlation equations for turbulent 
reacting flows. 

The construction of the "typical eddy" model is an attempt 
to utilize all the available information from a second-order 
closure analysis to define the joint pdf of the scalars. The 
degrees of freedom in the model, i.e., the number of delta 
functions, their strengths and their location in the phase space, 
are limited by the number of independent parameters available 
in a second-order closure calculation. The location of some of 
the delta functions is decided on physical grounds. The posi- 
tions of the other delta functions and the strengths of all the 
;delfa "functions at every point in the flow are determined from 
the values of the means and second-order correlations obtained 
from the solution of the transport equations. The method of 
construction of the delta function pdf model guarantees that the 
model reproduces the first-order and second-order moments of the 
flow. 


The complete second-order closure program for a three- 
species flow system provides information on 13 independent first- 
and second-order correlations - a , B , a ' 3 ' , a r y r , 3 * Y ' 
a’p' , 3 'P T , p'p' , h , h'h' , a ' h ' , B'h' , and p ' h' ; so 
the complete model can have no more than 13 parameters. Figure 1 
shows the model that has been selected after extensive studies. 

The model has exactly 13 parameters to be determined. The typical 
eddy is considered to have the following structure. For a frac- 
tion , the eddy contains only the species a . The cell has 

a sensible enthalpy h a and, since for a single species flow 
from a constant ent halpy' source the model must allow enthalpy 
fluctuations, h'h' , due to boundary conditions and/or viscous 
heating, a Ah has been added for half of the cell. For a 
single species flow, the enthalpy structure has just two para- 
meters (h a and Ah) to be calculated from the values of E and 
h'h' obtained from the transport equations. In the model cur- 
rently being used in the reacting shear layer programs, the same 
Ah is used for all the cells. For a, fraction £2 > the eddy 

contains only species 3 with an enthalpy h R , and for a frac- 

tion e^ , the eddy has only the pure species^ y and has the 
corresponding enthalpy hy . For fractions ejj , £5 , and eg , 

the eddy is assumed to contain a and 3 , a and y , and 3 

and y , respectively, in a state .of/ molecular mixedness. There 
appear to be two equally reasonable choices/for the proportions 
of the species in these cells. One can consider the species to 
be present in either equal amounts or. amounts proportional to ~ '3 
their average values at that point in the flow. The current 
program retains both of these options at the present time, and 
further test runs will indicate the more desirable model. Cor- 
responding average values are used for the enthalpies of the cells. 
Finally, an £7 fraction of the eddy contains all three species 
a , 3 , and y in a molecular mixture of proportions <1 5 <2 > 

and 1 - Kq - <2 • The enthalpy of cell 7 is hy . The model 
has 13 unknowns (e^ through eg , xq , <2 5 h a » h 3 5 h y » h 7 » 
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Figure 1. Complete three-species "typical eddy" model. 


1 4 










and Ah) note 


to be determined from the 13 avail- 


(l- 7 £i = 0 

able.. independent correlations by matching the moments of the 
model to the values predicted by the equations. The procedure 
for the construction of the "typical eddy" model is described 
in Appendix D, and further details can be seen in References 23 
and 24. The model can be constructed at any desired time from 
the second-order closure turbulence calculation in progress. 

Once the species and enthalpy distribution functions have been 
established, the corresponding distributions for the other 
scalars - temperature, density, and the reaction rate, k - can 
be constructed. The complete nonlinear Arrhenius rate expression 
can be used for constructing the k distribution, function; there 
is no need to expand the exponential term. All scalar correla- 
tions required to close the system of second-order correlation 
equations can now be calcualted from this joint probability 
density function. 

The construction of the complete "typical eddy" model 
described above requires the solution of a set of nonlinear 
algebraic equations. The nonlinear equations arise due to the 
density correlations and a simplified "typical eddy" model was 
proposed which was easier to implement for the initial testing 
of the concept of the model. 

Simplified "typical eddy" model .- The simplified "typical 
eddy" model does not consider the density fluctuation correla- 
tions in the construction of the model. The neglect_of the 

dens ity corr elat ion results in o nly nine moments - a , 3 , 

a'B’ , a'y' , B'y' , h , h f h f , a'h' , and B'h' - being 
available to formulate the model and, therefore, the simplified 
model can only have nine parameters. The model construction is 
briefly discussed in Appendix D. Further details of the construe 
tion and testing of the simplified "typical eddy" model are 
described in References 23, 24, and 25. The simplified model 
results in a linear set of equations for the determination of 
the cell sizes in the species probability density function, and 
analytical closed form solutions for the model have been obtained 
In contrast, one has had to resort to numerical solutions in 
the case of the complete three species model, although, some 
progress has been recently made towards obtaining analytical 
solutions of the nonlinear equations. This is discussed later 
in this section. 

The simplified model involves the neglect of the density 
correlations and, therefore, should be valid for flows where 
density changes are small. The simplified "typical eddy" model 
has been shown (ref. 26) to be satisfactory for low heat release 
reacting flows where the density changes could be neglected. 
However, just how far one could use this simplified scheme could 
only be tested by actual computations and comparison to experi- 
mental measurements. Our studies on chemical laser flows and 
hydrogen-air flames where significant density changes are in- 
volved suggest that the more complete "typical eddy" model has 
to be used, as was expected, for these flows. 
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The current 'reacting shear layer (RSL) computer program 
is designed to use the simplified "typical eddy" model or an 
even simpler model designated as the "secondary" model. The 
model used in a particular test run is conveniently selected 
as an input option. 

Secondary model .- In this model all third-order and higher- 
order scalar correlations are set to zero. This is designated 
as the "secondary" model, as distinguished from the primary 
model which is the "typ ical eddy" model. Second-order correla- 
tions such as T ! a' , k'a' , and others for which we do not 
carry transport equations, are now obtained up to second-order 
accuracy by expansion of explicit thermodynamic relations between 
the instantaneous variables. When the "typical eddy" model is 
used, these conversion relations are t not necessary as 

the pdf for p,T and k can be constructed directly from 

the pdf of and h using the thermodynamic relationships. 

The "secondary" model was initially used during the model develop- 
ment stage but some recent calculations (ref. 26) have indicated 
that, at least for some simple reacting flows, the two models 
may predict quite similar results. It can be shown quite easily 
that the "secondary" model is not correct in a number of limiting 
cases. For example, in the reaction end limit for a system of 
three species undergoing the one step reaction a + 3 y , 


a3 = 


a 3 + a ' 3 = 


a3’ + a'3 + a'3' - a'3' = 


a a'3 T + a'a' 3 + a'a'3' = 


a'a'3 1 = 


, / 0, in general 

Thus, setting the third-order scalar correlations to zero 
is not a proper model and the good agreement obtained between 
the results from the "secondary" model and the simplified 
"typical eddy" model must be simply coincidental. However, it 
is still an interesting, question why the results for some of the 
flows that we have studied are insensitive to the choice of 
the model for the higher-order scalar correlations, and work 
is continuing to explain these observations. 

The exothermic reacting flow calculations reported here 
have been carried out using the "secondary" model. These calcula- 
tions should be repeated later with the complete "typical eddy" 
model._ 


0 

0 

0 

0 


-a a'3' - 3 a' a * 
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Analytical solution for the complete "typical eddy" model . - 
Recent studies on chemical laser flows and hydrogen-air dif- 
fusion flames have shown that in flows involving large density 
changes the simplified "typical eddy" model leads to large errors 
and the more complete model including the density moments has to 
be used. We have completed the formulation of the nonlinear 
equations for the complete model and it appears likely that 
analytical solutions of these nonlinear equations can be obtained. 
Analytical closed form solutions for a "typical eddy" model for 
a two species mixing problem have already been obtained. The 
solution is very promising for it proves that for any statistically 
consistent set of moments, one can find a rational solution for 
the model parameters, that is, the strengths of the delta func- 
tions are positive and their location is inside the physically 
realistic region of the scalar phase space. Efforts are under- 
way to obtain similar analytical solutions for the three species 
model for the reacting flow problems. 

The reacting shear layer (RSL) computer program for second- 
order closure analysis of turbulent reacting flows will be able 
to incorporate the complete "typical eddy" model with very 
minor changes. 


Multi-Step Chemistry Procedure 


The reacting shear layer program and the "typical eddy" 
model were originally developed for use with three species 
(a , 3 , and y ) flowfields, and only permitted the single- 
step reaction a + -3 y . The concept of such an overall 
reaction may be entirely sufficient in some combustion problems, 
but for most problems it is necessary to use a set of elementary 
reaction steps involving many species to correctly describe the 
system and to obtain results in better agreement with experi- 
ments. This conclusion is particularly valid with respect to 
trace species. Further, the elementary reactions have been 
studied extensively and the reaction rate parameters are better 
known. Ideally, one would like to keep track of each species 
and the correlation- between each species and the other species 
as well as the correlation with certain flow properties. Un- 
fortunately, this desire introduces an enormous number of equa- 
tions which must be solved in an implicit fashion. The resulting 
computer time would be prohibitive to the point that the code 
would simply not be used for engineering calculations. There- 
fore, a system has been constructed for handling a finite num- 
ber of species and reactions within the present framework. 
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To understand the basic philosophy of this reaction scheme, 
it is useful to consider the expression for the rate of change 
of a species, a , due to chemistry in a situation in which 
turbulence may be present. It can be shown that the expression 
reduces to: 

+ aFnr*" + gk'a T + k'a' 3 * ] 

Turbulent 
Contribution 


= -kag - [ka ' 6 ' 

Laminar 

Contribution 


In the absence of turbulence, the rate of change of the 
species is simply given by the laminar expression. When turbulence 
is present, all of the other terms must be considered. The com- 
putation of these terms, and the other necessary correlations, 
is the reason the computer time becomes excessive. In the laminar 
case, such correlations are not required. Thus, the extension 
of laminar flow situations to hundreds of chemical species and 
reactions is straightforward since it only requires the addition 
of another equation for the conservation of mass of each added 
species. In the turbulent case, the addition of a new species 
requires the addition of at least three equations for the cor- 
relations with the primary flow variables plus an equation for 
the correlation with each of the other species in the flow. From 
a practical viewpoint, and the desire to maintain reasonable 
computation times, it was decided to "track" only three species: 
a , 6 , and y . 

In any diffusive reacting flow, one usually has two initially 
separate streams of fuel (a) and oxidizer ( 6 ) . These streams 
are now considered to be known mixtures of a number of chemical 
species. The streams mix and react following a large number of 
elementary reactions and form . a mixture of product species 
designated as y . It is assumed that (1) the composition of mix- 
tures a and 6 is fixed for the entire flowfield. There are 
no internal reactions among the species in these two reactant 
mixtures. (2) All the product species that compose y are 
molecularly mixed. There are internal reactions within the 
y species mixture. The composition of the mixture y varies 
at different points across the flowfield. 

Consider the hydrogen-air diffusion flame as an example of 
the multi-step chemistry procedure. To illustrate the procedure 
we consider only the two reactions given below. The complete 
set of reactions that are used in our studies of the hydrogen- 
air system are presented in Table 1. 

H 2 + OH ^=1 H 2 0 + H 

H 2 + 0 7 -^ OH + H 
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Let the a stream be pure hydrogen, H 2 ( a ) . The 3 

stream is a mixture of oxygen, 02 ( 3 ) and nitrogen, ^(3) 
of known, fixed proportions. The y product species contains 
H 2 , OH , H , 0 , O 2 , H 2 O , N 2 , etc. On the basis of 

the two assumptions mentioned earlier, the interaction between 
turbulence and ehcmistry only has to be taken into account for 
reactions between a and 3 , a and y , and 6 and y . 
Species within y are assumed to be molecularly mixed and, 
therefore, species correlations do not have to be considered 
for the internal reactions. For convenience, the reaction 
system given above is rewritten in the following form, tagging 
each chemical species with the mixture that it is a part of: 


H 2 (ci) 

h 2 (cO 

h 2 (y) 

H 2 (y) 


+ 


+ 


+ 


+ 


QH(y) 

0(y ) 

OH(y ) 

0 (y ) 








H 2 0(y) + H(y) 
OH(y) + H(y) 

H 2 0(y) + H(y) 
'OH(y) + H(y) 


one way reaction 


one way reaction 


The chemical source terms for the component species can 
now be written: 


dHp (a) 

— — = - k fl [H^Ta ) OH ( y ) + H*(cOOH» (y ) ] 

>■', - r k f2 [HjTa)OTyT + H^(a)0' (y)] 


dH^TyT 

dt 


kfi H^Ty)OHly) + k rl H 2 0(y)H(7y 
k f2 H 2 (y)0(y ) + k r2 OH(y)H(y) 

etc . 
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Transport equations are solved for the second-order correla- 


tions a' 3' , a’y' , and g'y* . The additional important 

assumption is made that quantities like HJ, (q ) OH ' ( y ) can be 
simplyy calculated from these correlations as being proportional 
to the local mean composition of the mixture. 

For example, H^a)OH'( Y ) = ^ ( a ° 'T'otffi. ) ^ 

l a • y J 

The validity of this assumption has to be examined further, but 
it should not lead to large error. 'The correct procedure, of 
course, is to solve separate equations for each such correlation, 
but this is impractical at the present time. 

The above procedure operates within the framework of three 
overall species mixtures, and all the models and equations 
developed earlier can be used. The only addition to the program 
is the solution of mean species conservation equations for all 
the elemental species. 


Numerical Solution of the Equations 


The use of the models for various third- and higher-order 
correlations, as described earlier, in the set of equations for. 
the means and second-order correlations results in a closed set 
of equations. Equations for the mean variables p , iL , E , 

a , and g and the second-order correlations u|u \ , h'h* , 

u|h' , a'g' , a'y' , g'y’ , u|a' , u|g* , u|p' , h' a ' , and 

h'g 1 are solved in the present program. The use of the basic 
shear layer assumptions leads to a set of 23 independent coupled 
parabolic partial differential equations. 

The computer program actually solves a total- of 29 equations 
for a three species system. The additional equations are for 

y , u’ y ' , v 1 y ’ , h' y 1 , a ' a ' and g ? g ' • The redundant equa- 

tions provide a check on the mass conservation in the program 
and on the accuracy of the numerical scheme. The numerical ^ 
integration of the equations is performed by a forward-time , 
upwind, finite-differencing scheme. The nonlinear terms are 
handled by quasi-linearization, that is, by evaluating a portion 
of these terms at the known position, leaving only a linear term 
containing one of the' unknowns,-. The linearized finite-difference 
equations are solved in implicit fashion using the general tri- 
diagonal algorithm. However, instead of solving a single large 
matrix, which would be very time-consuming, the equations are 
grouped into smaller matrices and the system is solved in 10 
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separate passes. The equations solved in the various passes are 
detailed in the table below. The variables are updated only 
at the end of all the 10 passes. For multi-step chemistry prob- 
lems, the first pass solves the conservation equations for all 
the elemental species. 


Table RSL Passes 


Pass No. Unknowns 



a , 3 , 





1 

y or 

a . 

l 


2 

u , <u r v 

' > 




3 

-C 

> 

V 

r* 

ix: 

’> , < 

P 

' V’ 

' > 

4 

<u'u f > , 

<v' V* 

> 

3 

<w'w r > 

5 

<h'h’> , 

<u'h' 

> 

3 

<p'u' > 

6 

<u ' a ' > , 

<u' 3’ 

> 

3 

<u ' y 1 > 

7 

A 

< 

P 

V 

v* 

<v' 3’ 

> 

3 

A 

>- 

> 

V 

8 

r» 

A 

a 

V 

<h ' 3 ' 

> 

3 

A 

>- 

xi 

V 

9 

<a'a'> , 

<a ' 3 ' 

> 

3 

A 

>- 

3 

V 

10 

<3 ’ 3 ’ > , 

<3 ' y ' 

> 



+ continuity 

equation for v 






The program has the capability of handling both fixed and 
free shear, layer flows for planar and axisymmetric geometries. 

The' initial .profiTes-can -be. provided using a card or tape input 
or the program can -generate appropriate smooth initial profiles 
given the properties of the two streams. ' If the turbulence 
properties at the initial station are not known, a "spot” of 
turbulence is input at the initial station to start the tur- 
bulent calculations. The turbulence profile is a smooth bell- 
shaped profile of specified small amplitude and width. It has 
been shown that the calculated results after approximately 10 
diameters downstream are insensitive to the amplitude and width 
of the initial turbulence profile. The program uses a fourth- 
order polynomial expression to calculate the temperature dependence 
of the specific heat of various species. The molecular trans- 
port properties y , y* , k , and D are calculated using a 
power law expression. 

Additional details of the program operation for a number 
of the test runs are discussed later in the report. Appendix C 
contains a description of the RSL program input deck and the 
details on the use of the various available options. 
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MODEL AND PROGRAM VERIFICATION - NONREACTING FLOWS 


Introduction 


The principal objective of developing a higher-order closure 
turbulence model than the first-order closure eddy-viscosity 
models is to obtain a more "universal" model; a model which can 
be used for analysis of data for a wide class of flowfields with 
the use of the same invariant set of model constants. Such a 
model can then be used with a greater degree of confidence for 
predictive calculations of flow problems for which experimental 
data are not available or are difficult to obtain. Second-order 
closure models have shown significant promise of being such ; 'a 
predictive tool. 

The fluid mechanical turbulence models being used in the 
A.R.A.P. second-order closure program have already been extensively 
tested by comparison of program predictions to experimental data 
in a variety of basic flow geometries. Previous publications 
have discussed studies on flat plate boundary layers (refs. 15,22), 
the planetary boundary layer (ref. 16), and two-dimensional wake's , 
aridjjets (ref. 26). The two-dimensional reacting shear layer' 

(RSL) code was extended to axisymmetric flows under the NASA . .. 
contract, and model and program verification studies for axi T 
symmetric jets and wake flowfields and the mixing of two dif--; 
ferent species in a shear layer geometry have been carried-'out 
in the following sample cases. . ■ 

Axisymmetric Free Jet 


'The RSL program predictions for an axisymmetric free jet 
were compared to the experimental measurements of Wygnanski and 
Fiedler (ref. 27). Detailed measurements of the mean axial 
velocity profile and the Reynolds stresses were made in the self 
similar region of the flowfield. The measurements are designated 
as Test Case 18 in data compiled in the 1972 NASA Free Turbulent 
Shear Flows Conference (ref. 28). 


The A.R.A.P. models have already been tested against this 
flowfield in previous simpler nonreacting programs (refs. 17,29) 
and the present repetition of the test case provides a good 
consistency check for the RSL code. The. computations are started 
using a simple linear (arbitrary) profile for the mean velocity 
and a triangular distribution for the Reynolds stresses with 

u’ u"’ = v* v ' = w * w 1 = 2 | u’ v ’ j at the initial station. The 
turbulence scale A is determined from the turbulence kinetic 
energy q2 profile. The model constants for the diffusion 
terms are P 2 = -0.3 and V c = + .3 . The constants P-^ and P^ 

are set at zero. 
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Comparisons between RSL program model predictions and ex- 
perimental data are shown in Figures 2, 3, and 4. Figure 2 
compares the mean axial velocity profiles and indicates ex- 
cellent agreement. Figure s 3 and 4 s how the results for the 
normal Reynolds stresses u ' u * and v'v' . The predicted 
values are somewhat low on the centerline. The agreement can be 
improved for this flow by incorporating the mean strain terms 
in the pressure correlation models but the values of the cor- 
responding model parameters have not yet been evaluated. Further, 
to obtain a model that is reasonably good for a wide range of 
flows, it is best not to fine-tune it to a particular flow. 

The parameters for the complete pressure models will be selected 
to improve the model predictions for both incompressible and com- 
pressible flows, and then these basic test flows will be re- 
examined. The RSL program results are consistent with previous 
A.R.A.P. studies. » 


Axisymmetric Wake 


The second axisymmetric test flow was selected to be the 
axisymmetric wake measurements of Chevray (ref. 30). A.R.A.P. 
model predictions have also been compared against other axi- 
symmetric wake measurements in Reference 29. Chevray 
measured the mean velocity u and the Reynolds stress u' v* 
at several axial stations in the wake behind a six-to-one prolate 
spheroid at a Reynolds number, based on free stream velocity 
and body length, of 2.75 x 10°. The data is NASA Test Case 15 
in Reference 28. The measurements at the axial station three 
diameters downstream of the trailing edge of the body were chosen 
as the initial conditions for the model calculations. This was 
assumed to be far enough downstream that our assumption of no 
axial pressure gradient in the calculations would be valid. The 
initial conditions for the program run used the measured data 

for u and u 1 v' at = 3 . The normal stresses at this 

initial station were not known and were assumed to be 

u' u' = v’v* = w ’ w T = 2 | u ' v ’ | . The scale A is :c„ohstant across 

the flow and is determined from the turbulence energy q 2 • n 
profile. The model constants are again P_ = -0.3 , V = +0.3 
and Pq = P3 = 0 • d C 

Figure 5 shows the decay of the maximum velocity defect 
Wp . The model predictions are reasonably good, though the data 
appears to decay somewhat faster than the analytical predictions. 

The maximum shear stress u'v^ is plotted in Figure 6. The 

model predicts a higher initial buildup of the shear stress 
before it begins to decay and is then in good agreement with the 
data. The predicted profiles of the mean velocity and the shear 

stress at axial stations ^ = 9 and 18 are compared to the measure 
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Figure 2. Mean axial velocity profile in a free jet in the 
self-similar region. 
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A.R.A.P. prediction 



Figure 3* Profile of the axial velocity fluctuations in a 
free jet in the similarity region. 
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A.R.A.P. prediction 


O Wygnanski & Fiedler 27 
data 



Figure 4 . Profile of the radial velocity fluctuations in a 
free jet in the similarity region. 
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x/D 

Figure 5. Variation of the mean velocity defect with distance 
downstream of the body. 



ments in Figures 7 and 8, and show good agreement. Previously 
reported calculations of this flowfield using a dynamic scale 
equation (ref. 17) have shown somewhat better agreement with 
the data. The RSL calculations using a constant scale A are 
again consistent with previous (ref. 29) A.R.A.P. calculations 
using an incompressible program, and verify the operation of the 
reacting shear layer (RSL) program. 

The RSL program has been similarly checked for other flow 
configurations in previous publications, for example, flat plate 
boundary layers (ref. 22), two-dimensional jets and wakes 
(ref. 26), etc. 


Uniform Density Shear Layer Mixing of Two Species 


Recent measurements in a two-species mixing layer have 
been reported by Konrad (ref. 31): These measurements are part 

of a very elaborate and detailed study of this basic flowfield 
by the group under the direction of Professor Roshko at Caltech. 
Konrad has made measurements of the mixedness correlation 

a'g'/aB or a’a'/a(l-a) in a shear layer of velocity ratio 
= 0.38 and consisting of two streams of different species but 
the same density. Figure 9 shows a comparison of the RSL 
program predictions for the flowfield with the experimental 
measurements. There is good agreement in the central turbulent 
region of the flowfield. The theoretical predictions for the 
edge regions of the flow are not symmetrical, which is somewhat 
surprising and is believed to be due to the difficulty of a 
accurately calculating the term a’g'/ag when a or g is 
very small. Further, the calculations were carried out using 
the secondary model and the neglect of the third-order density 
correlations may not be justified in this flow. The calculations 
should be repeated with the complete "typical eddy" model. The 
disagreement between the theoretical predictions and the experi- 
ments may also be due to the intermittent nature of the flow and 
the observed large "coherent" structures. This is the only 
available data on this flow configuration with detailed turbulence 
measurements, and the comparison to the theoretical predictions 
to this flow have been simply included in this report' as an area 
that requires further study. The theoretical models do not in- 
clude intermittency effects at the present time. 
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Konrad 

Experimental data 



Figure 9. Profile of the unmixedness correlation - in a shear layer 






Conclusions 


The most important feature that must be pointed out is 
that all these computations and others that have been referenced, 
have been carried out with basically the same set of model 
parameters. The results demonstrate that our second-order 
closure techniques can predict the mean flow and details of the 
turbulence correlations for a wide class of nonreacting flows 
and provide a. level of confidence in the use of, these models 
and the programs in a predictive capacity for other flows. The 
details of the specific agreement with one particular experiment 
are not as important in this effort as the approach to the 
development of a universally valid model. The reacting shear 
layer (RSL) program is, therefore, at least correct for these 
limiting cases and now provides the capability to study multi- 
species mixing and reacting flows. Tests now have to be carried 
out on these flows to verify the rest of the program and to 
verify the models developed for the scalar correlations. This 
effort is now underway. Computations are being carried out for 
two species shear layer mixing (the initial preliminary results 
were shown previously in this section), hydrogen-air diffusion 
flames (see next section), DP chemical lasers (ref. 10;), and 
erosive burning of solid propellants (ref. 32). However, this 
work is just beginning and many tests and comparisons remain 
to be made to thoroughly check the models and the programs. 


STUDIES ON HYDROGEN-AIR DIFFUSION FLAMES 


There are at the present time very few detailed measurements 
in turbulent reacting flows which can be used to test the assump- 
tions and the models in a complete second-order closure analysis. 
Measurements are required for various second-order correlations 
and especially the scalar correlations and the scalar probability 
density functions to guide the development of these more sophis- 
ticated turbulence models for reacting flows. These experiments 
are quite difficult, but many measurements of this kind are now 
being planned by various groups and some results should become 
available in the near future. At present, the only systematic 
set of data for testing theoretical models are the detailed measure 
ments of Kent and Bilger (ref. 11 ) on laboratory hydrogen-air 
diffusion flames. The A.R.A.P. reacting shear layer (RSL) program 
has been used to study this flowfield. 

A schematic of the Kent-Bilger experimental setup is shown 
in Figure 10. The measurements of temperature, velocity, and 
species compositions were carried out in a coaxial mixing system 
consisting of a 7.62 mm diameter hydrogen-jet exhausting into an 
air stream flowing through a 305 mm square cross-sectional tunnel. 
The flame was initiated by a spark. Measurements were reported 
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Length = 1.8 m 



Figure 10. Kent-Bilger experimental set-up. 


at nominal jet to external stream velocity ratios of 2, 5, 8, and 
10. The most extensive measurements were for a velocity ratio of 
10 and the theoretical calculations have been made for only this 
case. The actual velocity ratio for the measurements was 11.8. 

The measurements indicate the presence of a fairly strong axial 

. pre ssu ne — gradient . Rhode s e_t._._a.l_. (ref . 3_) hav^ e shown t: hp imnnr ’t- an.t 

effect of the pressure gradient on the mean velocity field. Al- 
though the RSL program has the capability of imposing an external 
pressure gradient (for example, on a boundary layer flow) at the 
present time, the program does not treat a confined flow in a con- 
sistent manner. This capability will be added in the near future. 
For the calculations reported herein, the axial pressure gradient 
was neglected. The effect of the inclusion or the neglect of the 
pressure gradient on the predictions of the temperature and the 
mean species profiles is not known. 


Measured velocity profiles at the nozzle exit plane were 
available, but to avoid the initial steep boundary layer profiles, 
the RSL calculations were started at x/D = 6 . The flow properties 
for all the mean variables at this axial station were obtained from 
a NASA Langley computer program which uses a two-equation turbulence 
model. This code generated the profiles given in Figure 11, which 
were the initial conditions for all the calculations discussed here. 
A narrow spot of turbulent velocity correlations 




was intro- 


duced at the initial station with amplitude of u'u r = v'v' = w'w' 


_2 

= 2u'v' = .01 Au . All the other second-order correlations were 
set zero at the initial station. 


A number of different program assumptions were investigated 
as part of the studies on the Kent-Bilger hydrogen-air diffusion 
flame . As most of the measurements were in the far downstream 
region of the flow, a simpler and faster "superequilibrium" procedure 
(ref. 18 ) was used to calculate the species correlations in a 
number of the test runs. In this procedure, the transport equations 
for the species correlations are not solved, but instead, the 
species correlations are calculated from algebraic relations ob- 
tained by balancing the production and dissipation terms in the 
complete transport equations. The convective and diffusive terms 
in the equations can be assumed to be small In the downstream region 
of the flow. The species correlations can then be written in 
terms of the mean flow variables, and have the form. 


a 1 ft 


= C ,A 2 
<P 


3a 

7y 


3 ft 

3y 


where ft may be any variable ( 0 , u, h, p, etc.) and C<j, is a con- 
stant that can be evaluated from the turbulence model parameters. 

In the studies discussed in this report the superequilibrium 
procedure was used only for the species correlations. Other 
turbulence correlations can also be evaluated with the use of the 
superequilibrium assumption. This aspect of the investigation 
should be studied in the future. The results of such superequili- 
brium calculations have to be compared to results obtained by use 
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Figure 11. Initial conditions for A.R.A.P. calculations 





of the transport equations for all the correlations which is the 
normal mode of program operation. Such a comparison was attempted 
here, but was not completed due to problems which are discussed 
later in this section. 


The turbulent macroscale A was calculated from the mean 
axial velocity profile for the hydrogen-air diffusion flame 
calculations. In a reacting flowfield, the profile of scalar 
variables, for example, the mean temperature is quite different 
from the mean velocity profile. Therefore, it is quite likely 
that the turbulence macroscales for different variables will also 
be different. However, there is noQ convenient procedure to 
accurately determine these different scales at the present time 
and all other analyses use the same scale for all the variables. 

At the present time A.R.A.P. is engaged in formulating a new and 
consistent procedure for determining these scales. In general, 
the scales for different variables will be different and will vary 
across the flowfield. 


The effect of the different scales was investigated in a very 
simple manner in the Kent-Bilger studies .^CTwo^I^ies_wereJ^qns*id^red. 
In one the scales J'or all the variables were the same and the scale 
-was. ...evaluated from ...the mean axial velocity profile. This scale 
is constant across the flowfield, but varies with the profile in 
the flow direction. In the second case, the scale for all the 
scalar variables is set equal to one-half the scale for all the 
turbulent velocity correlations. The 1/2 figure is an arbitrary 
estimate from the profile shapes at the initial station. The 1/2 
ratio is maintained for the entire flowfield and is likely to be 
in error, but a more complicated formulation did not appear 
justified at this time. The detailed scale equation formulation 
currently underway should lead to a better understanding of the 
physics of the problem and then a more appropriate selection of the 
scales may be possible. Thus, in the current studies we investi^?. 
gated. 


^scalar A u 
A scalar = 1//2 A u 


All the calculations used the secondary model for the scalar 
correlations. The calculations using the complete "t ypical eddy" 
model with the density correlations will be done 1 a t er_af ter the 
complete model is incorporated into the RSL program. 


The chemical rate constants used for the multi-step hydrogen- 
air chemistry are listed in Table I. These are in gerferal agree- 
ment with the rates used by other investigators (refs. 3,33)- 


There is one additional empirical parameter that has to be 
specified in the current reacting runs using the RSL program. 

This is the limiting value for the unmixedness correlation in the 
chemical source term. The mean chemical source term for the 
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elementary reaction a + 8 — >-Y can be written as. 


w 


— 1 r ^ /M 


a 


kp 


a8 1 


+ + 2 sigl +■ 


ap 


+ k T terms 


8jP_! + a’g* 

BP a8 


The total set of correlations inside the square brackets are 
termed the unmixedness correlation. For an infinitely fast reac- 
tion rate, the unmixedness correlation = -1. and the species 
a and 8 cannot coexist on the molecular scale. This is the 
limiting value for k = °° . However, for large k , the unmixed- 

ness correlation approaches -1 , but is not equal to it. In princi- 
ple, the unmixedness correlation for large k can be calculated 
by solving the transport equations for the individual correlations. 
However, this procedure is subject to errors when the reaction 
rate is fast and it becomes necessary to calculate the limiting 
value of the unmixedness correlation bv some different procedure. 

A fast chemistry limit has been derived for constant density 
reacting flows (see Appendix B) and the procedure is being extended 
to exothermic, variable density reacting flows. Alternately, the 
limiting value can be input as an empirical parameter. In the 
current set of calculations for the Kent-Bilger hydrogen-air diffu- 
sion flame studies, the alternate procedure of an empirical limiting 
value has been used. The unmixedness correlation is calculated 
at every point in ’/the flow using the transport equations or the 
super equilibrium procedure. If the value of the unmixedness is 
smaller than the limiting value, it is set equal to the limiting 
value at that point in the flowfield. In the hydrogen-air dif- 
fusion flame problem, the reaction rates are very fast and the 
limiting values of the unmixedness correlation were found to be of the 
order of -0.996 to -0.999. Thus, the net turbulent reaction rate 
is only a very small fraction (.001 to .004) of the reaction rate 
that would be predicted if the effects of the turbulence were 
neglected and the flow was assumed to be molecularly mixed. It 
is, therefore quite obvious that the turbulence-chemistry inter- 
action is of great importance here. However, due to the fast 
reaction rates of this problem, it is difficult to accurately com- 
putei'the unmixedness correlation from the basic transport equa- 
tions. The calculation of the unmixedness effects will be con- 
siderably simpler for slower reaction systems, as the need for 
extreme accuracy in the numerical solution of the transport equa- 
tions for the correlations will be obviated. 


The results obtained for 


A = A- are considered first, 

scalar u 


Figures 12 and 13 compare the results for the temperature and 
species profiles for the RSL predictions with the Kent-Bilger 
measurements at x/D = 40 . The calculations use t.he superequili- 
brium procedure and use an unmixedness limit of -0.9995. The 
peak temperatures are in reasonably good agreement but the calculated 
temperature profile is significantly wider than the measured profile 


38 





oOA Experimental data 

Turbulent chemistry 

calculations 



Superequilibrium procedure. 


at this position. The species profiles also show that the flame 
location is predicted at a somewhat larger radius, than in the ex^ 
periments. The shape of the oxygen profile is steeper in the 
measurements but the hydrogen profile is not too -different if one 
compensates for the difference in flame location. The measured 
flame thickness (region of overlap of hydrogen and oxygen) is 
larger than the predictions which suggests. that the unmixedness 
limit should be set even closer to -1 . Computations made with the 
limit set at -.998 and -.999 predicted even smaller reaction zone 
thickness. The results at x/D = 80 are shown in Figures 14 and 
15 and the same comments apply although the profile, shapes for the 
temperature and the oxygen are now in better agreement with the 
data. 


To show the importance of including the turbulence-chemistry 
interaction effects, a computation was attempted using a ’’laminar 
chemistry" formulation. In this case the unmixedness correlation 
is completely neglected in the calculation of the reaction source 
terms and the flow is assumed to be perfectly mixed on the molecular 

scale-Ao Then, w = kp 2 ag , and is analogous to the source term, 
for laminar flow a which is the reason for the name used for this 
approach. Unfortunately, the run gets into severe ^problems with 
the prediction of very large values of h ’h ’ and consequent 
errors in the calculation of the mean temperature and density, 
and the comparison could not be made for A scalar = A 5 case. 

( ] > \ r\ 

A' test calculation which,does not use the superequilibrium 
assumption for the species correlations was also attempted, but 
was not successful due to problems in the multi-step chemistry 
procedure. The forward step size in the computation becomes 
very small due to backups in the minor species diffusion and 
chemistry calculations and the computation could not be completed. 
This problem can be corrected but up to now we have not been able 
to devote the necessary effort. 

The problems experienced in the "laminar chemistry" run for 

A , = A- (too large h ’h ’ ) and the predicted large profile 

scalar u 

widths for the temperature and the species in the "turbulent 
chemistry" runs suggested that the macroscale for the scalar 
quantities probably should be smaller. As discussed before, on 
the basis of the initial profiles of the variables, we selected 
Ascaiar = A- make a sep i es test computations for the 

Kent-Bilger diffusion flame studies. 

Figures 16 and 17 compare the program predictions to the 
measurements at x/D = 40 for this scale assumption. The calcula- 
tions uset the superequilibrium assumption . ; The unmixedness cor- 
relation limit is set at -.9-9.63 after a number of tests. The 
results are quite sensitive to the choice of the limit as is shown 
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later. Figure 16 shows the results for the mean species with the 
use of both the "laminar chemistry" (.perfect mixing.) and the 
"turbulent chemistry" formulations. It is quite, clear that there 
is much better agreement with the experimental data with the -in- 
clusion of the scalar c orr e lat ions'jin the reaction, source terms. 
When the species and density fluctuations are neglected in the 
"laminar chemistry" approximation, the flame thickness is very 
small. The turbulence-chemistry interaction effects have to be 
included to correctly predict the observed turbulent flame thick- 
ness. /However, there still is significant disagreement between 
the computer predictions and. the data. The predicted flame 
position is at a larger radius than the measurements indicate and 
the overlap between the hydrogen and oxygen profiles is now larger 
than in the experimental observations. These differences require 
further testing of the models and the assumptions used in these 
calculations, for example, the neglect of the axial pressure 
gradient, the macroscale assumption, the superequilibrium assump- 
tion, etc. 

Figure 17 shows the results for the temperature profile. With 
the use of the "laminar chemistry" approximation the peak tempera- 
ture is high. The use of the "turbulent chemistry" formulation 
predicts temperatures that are much more consistent with the ex- 
perimental values, although the predictions are somewhat lower. 
Figures;)l6 and 17 clearly demonstrate that for calculations using 
finite rate chemistry, the significant features of turbulent dif- 
fusion flames can only be predicted when the effect of scalar 
fluctuations is included in the reaction source term computation. 
The figures also show that the RSL program predictions are in 
fairly good agreement with the Kent-Bilger hydrogen-air diffusion 
flame data at x/D = 40 when all the previously discussed assump- 
tions and empirical constants are used. However, when the.- predic- 
tions are compared to the data at x/D = 80 , the agreement is not 
satisfactory. Figures 18; and 19 show the results at this axial 
position. The significant disagreement is in the reaction zone 
thickness and in the peak temperature. It appears that the un- 
mixedness limit should be closer to -1 at this station than at 
x/D = 40 . This is not unreasonable, in fact, our studies on a 
simple mixing layer with chemical reaction (ref. 26) show that 
the unmixedness correlation asymptotically approaches -1 as x 
increases. However, it is difficult to incorporate this variation 
of the unmixedness limit in the theory at this time without ad- 
ditional empiricism. 

Figure 20 shows the sensitivity of the predictions to the 
choice of the unmixedness limit . The calculations are compared 
at x/D = 40. The calculations with an unmixedness limit of 
-0.996 are in best agreement with the' measurements at this axial 
position. It is seen that the unmixedness correlation has to be 
calculated to an accuracy of the order of 10“3 . The use of a 
fast chemistry limit may prove to be the most accurate procedure 
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Figure 18. Mean species radial profile in hydrogen-air diffusion flame. Secondary model. 
Superequilibrium procedure. 
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Figure 20. Sensitivity of mean temperature profiles to unmixedness limit. Secondary model 
Superequilibrium procedure. 




for calculating the unmixedness correlation limit for the fast 
hydrogen-air reactions. 

In summary, the results of these calculations have demonstrated 
the importance of including the species and temperature fluctua- 
tion effects in calculations of turbulent diffusion flames with 
finite rate chemistry. These effects are included in a second- 
order closure analysis of the flow-field as is being developed by 
A.R.A.P. However, a considerable amount of work still remains to 
be done before one can obtain definitive comparisons between the 
reacting shear layer (RSL) program predictions and experimental 
data such as the measurements of Kent and Bilger. The importance 
of the selection of the correct macroscale for the scalar cor- 
relations and that for the turbulent velocity correlations has 
been demonstrated. The development of a scale equation formulation 
based on the two-point correlation function should be helpful in 
determining the proper scales. The empirical choice of the un- 
mixedness limit may be avoided by the development of suitable 
fast chemistry limits, or it may become unnecessary when the trans- 
port equations are accurately solved with the complete "typical 
eddy" model. The calculations reported here with the secondary 
model and the superequilibrium procedure have to be repeated with 
the use of the full solution of all the transport equations and 
with the "typical eddy" model for the scalar correlations. The 
concept of the simplified "typical eddy" model has been success- 
fully verified in low heat release reaction flows. However, the 
use of the simplified "typical eddy" model was not successful in 
the diffusion flame studies as substantial errors in the mean 
density and the density correlations occur due to their neglect 
in the model construction. Significant progress has been made 
recently toward the development of the complete model which in- 
cludes the density fluctuation correlations. 

SUMMARY AND CONCLUSIONS 


A complete second-order closure program for the study of 
turbulent reacting shear flows has been developed. The program 
is capable of handling multi-species, multi-step chemical reac- 
tions and includes turbulence-chemistry interaction effects in the 
analysis. Calculations for a hydrogen-air axisymmetric diffusion 
flame show the importance of including the scalar correlations of 
species and temperature in calculating the chemical source terms 
for finite rat e chemis try calculations. It is_n ecessar.v to include 
the turbulence-chemistry interaction to correctly predict important 
features of turbulent diffusion flames such as the thick reaction 
zone . 


The current studies used a number of assumptions and simplifi- 
cations, such as the use of the superequilibrium formulation, the 
use of the secondary model for the scalar correlations, the arbitrary 
specification of the turbulent scale for the scalars, and the 
unmixedness limit, etc. These approximations became necessary due 
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to a variety of reasons. The superequilibrium assumption was 
used simply because of numerical problems in the multi-step 
chemistry code and it is probably the simplest to eliminate. 

The secondary model had to be used because of the inadequacy of 
the simplified "typical eddy" model to handle variable density 
flows and this assumption can only be corrected when the complete 
model is developed. The calculations led to the discovery that 
it may be necessary to specify different macroscales for scalar 
variables in reacting flows. This has to be examined in greater 
detail before a final conclusion can be reached. The possibility 
was anticipated in the very early model development work (ref. 18) 
but was not necessary in any of our previous studies. The unmixed- 
ness limit can now be calculated using the procedure derived in 
Appendix B. The current studies should, therefore, be considered 
as preliminary and have served to define a number of problem areas 
that need further investigation. The development of the "typical 
eddy" model and the questions regarding the turbulent scale 
specification are probably the two most important areas that will 
change the results of the calculations. 

The A.R.A.P. reacting shear layer (RSL) program developed 
under this research effort is an operational program and will 
provide the basic capability to carry out second-order closure 
calculation of turbulent reacting flows. The new models developed 
in the future can be easily incorporated into the program. Our 
studies have already demonstrated some of the important effects 
in nonequilibrium turbulent flowfields that can only be properly 
analyzed with the use of a second-order closure program. 


REFERENCES 


1. Donaldson, C. duP.; and Hilst, G.-R.: Chemical Reactions in 

Inhomogeneous Mixtures: The Effect of the Scale of Turbulent 
Mixing. Proceedings of the 1972 Heat Transfer and Fluid 
Mechanics Institute, Stanford University Press, Stanford, CA, 

pp. 256-261. 

2. Spalding, D. B. : Mixing and Chemical Reaction in Steady Con- 

fined Turbulent Flames, 13th Symposium (International) on 
Combustion, the Combustion Institute, Pittsburgh, PA, 1971, 
pp . 6 4 5—6 57 - 

3. Rhodes, R. P.; Harsha, P. T.; and Peters, C. E. : Turbulent 

Kinetic Energy Analyses of Hydrogen-Air Diffusion Flames, 

Acta Astronautica, Vol. 1, 1974, pp. 443-470. 

4. Bray, K. N. C.; and Moss, J. B.: A Unified Statistical Model 

of the Premixed Turbulent Flame. AASU Report No. 335, Univer- 
sity of Southampton, November 1974. 

5. Libby, P. A.: One Turbulent Flows with Fast Chemical Reactions. 

III. Two-Dimensional Mixing with Highly Dilute Reactants. 
Combustion Science and Technology (Special Issue on Turbulent 
Reactive Flows), Gordon and Breach Science Publishers, NY 

(in press ) . 


51 



6 . 


Lockwood, F. C.; and Naguib, A. S.: The predictions of the 

Fluctuations in the Properties of Free, Round-Jet, Turbulent, 
Diffusion Flames, Combustion and Flame, Vol. 24, 1975, pp . 109- 


7. Borghi, R.; and Moreau, P.: Turbulent Combustion in a Premixed 

Flow. Paper presented at Fifth International Colloquium on 
Gasdynamics of Explosions and Reactive Systems, Bourges, France 
September 1975- 

8. Spalding, D. B.: A General Theory of Turbulent Combustion, 

the LaGrangian Aspects. Paper presented at AIAA 15th Aero- 
space Sciences Meeting, Los Angeles, CA, AIAA 77-141, January 
24-26, 1977. 

9. ' Kewley, D. J.: A Model of the Supersonic HF Chemical Mixing 

Laser Including Turbulence Effects on the Chemistry. Paper 
presented at D.G.L.R. Symposium on Gasdynamic and Chemical 
Lasers, Koln, W. Germany, October 1976. 

10. Varma, A. K.; Fishburne, E. S.; Beddini, R. A.: Modeling of 
i Turbulent Mixing and Reactions in Chemical Lasers. AFOSRi ■ 

l TR-77-0584, January 1977- 

11. /Kent, J. H.; and Bilger, R. W. : Turbulent Diffusion. Flames . 

{Fourteenth Symposium (International) on Combustion (The Com- 
bustion Institute, Pittsburgh, PA), 1973, pp. 615-6 25V ' ^ g 

12. Launder, B. E.; and Spalding, D. B.: Mathematical Models of... 

Turbulence. Academic Press, 1972. 

13. Bradshaw, P; : The Understanding and Prediction of Turbulent' . 

Flow, Aeronautical Journal, Vol. 16, 1972, pp. 403-418. 

14. Saffman, P. G.: A Model for Inhomogeneous Turbulent Flow, 

Proceedings of Royal Society of London, Vol. 317, 1970, 
pp. 417-433. 

15. Varma, A. K.; Beddini, R. A.; Sullivan, R. D.; and Donaldson, 

C. duP..: Application of an Invariant Second-Order Closure 

• Model to Compressible Turbulent Shear Layers. Paper presented 
at AIAA 7th Fluid and Plasma Dynamics Conference, Palo Altp, 

CA, June 17-19, 1974, AIAA 74-592. ; 

16. Lewellen, W. S.-; Teske, M. E.; and Donaldson, C. duP.: Z Tur- 

bulence Model of Diurnal Variations in the Planetary Boundary 
Layer. Proceedings of the 1974 Heat Transfer and Fluid/ 
Mechanics Institute, 1974, pp . 301-319.- />: 

/ ' 

17. Lewellen, W. S. : Use of Invariant Modeling. A.R.A.P.J .Report 

No. 243, April 1975. ./ i : * r f 

18. Donaldson, C. duP.: .. Atmospheric Turbulence" and • ttfe 'Dispersal 
of Atmospheric Pollutants. EPA-R4-73-0l6a, 1973. 


52 



19. Fishburne, E. S.; Beddini, R. A.; and Varma, A. K.: The Com- 

' putation of Afterburning Rocket Exhaust Plumes. A.R.A.P. 

Report No. 283, July 1976. 

20. Varma, A. K.: ‘Second-Order Closure of Turbulent Reacting 

Shear Flows. A.R.A.P. Report No. 235, February 1975. 

21. Lewellen, W. S.; and Teske, M. E.: Prediction of the Monin- 

Obukhov Similarity Functions from an Invariant Model of 
Turbulence, J. Atmos. Sci., Vol. 30, 1973, pp. 1340-1345. 

22. Varma, A. K.; Beddini, R. A.; Sullivan, R. D.; and Fishburne, 

E. S.: Turbulent Shear Flows in Laser Nozzels and Cavities. 
A.R.A.P. Report No. 231, October 1974. 

23. Donaldson, C. duP.: On the Modeling of the Scalar Correlations 

Necessary to Construct a Second-Order Closure Description of 
Turbulent Reacting Flows. Turbulent Mixing in Nonreactive 

and Reactive Flows, S.N.B. Murthy, ed.. Plenum Press, (New York) 

1975, pp. 131-162. 

24. Donaldson, C. duP.; and Varma, A. K.: Remarks on the Con- 

struction of a Second-Order Closure Description of Turbulent 
Reacting Flows, Combustion Science and Technology, Vol. 13, 

1976, pp. 55-78. 

25. Varma, A. K.; and Donaldson, C. duP.: Development of a Second- 

Order Closure Model for Computation of Turbulent Diffusion 
Flames. NASA CR-132524, June 19.74. 

26. Varma, A. K.; Beddini, R. A.; and Fishburne, E. S.: Second- 
Order Closure Analysis of Turbulent Reacting Flows. Proceedings 
of the 1976 Heat Transfer and Fluid Mechanics Institute, 

1976. 

27- Wygnanski, I.; and Fiedler, H. : Some Measurements in the Self- 

Preserving Jet, J. Fluid Mech., Vol. 38 , pt . 3, September 18, 
1969, PP. 577-612. 

28. Birch, S. F.; and Eggers, J. M. : Free Turbulent Shear Flows, 

Vol. I and II, Conference Proceedings, NASA SP-321, 1973. 

29. Lewellen, W. S.; Teske, M. E.; and Donaldson, Gu duP.: 
Application of Turbulence Model Equations to Axisymmetric 
Wakes, AIAA J., Vol. 12, No. 5, 1974, pp. 620-625. 

30. Chevray, R.: The Turbulent Wake of a Body of Revolution. 

Trans. ASME, Ser. D. : J. Basic Eng., Vol. 90, No. 2, June 1968, 

pp. 275-2 84. 


53 


31. Konrad, J. H.: An Experimental Investigation of Mixing in 

Two-Dimensional Turbulent Shear Plows with Applications to 
Diffusion-Limited Chemical Reactions. Ph.D. Thesis, California 
Institute of Technology, Pasadena, CA, 1976. 

32. Beddini, R. A.; Yarma, A. K.; and Fishburne, E. S.: A Pre- 

liminary Investigation of Velocity-Coupled ■ Erosive Burning, 
Proceedings of the 13th JANNAP Combustion Meeting, Vol. II, 
December 1976, pp. 385-^06. 

33* Pergament, H. S.; and Fishburne, E. S.: Influence of Buoyancy 

on Turbulent Hydrogen/Air Diffusion Flames. Paper presented 
at Fluid Mechanics of Combustion Processes Conference, Central 
States Section, The Combustion institute, March 28-30, 1977. 


5 \ 
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APPENDIX A. 


General Conservation Equations for a Turbulent 
Multicomponent Reacting System 


Continuity Equation 


p t + ( pu ) ^ = 0 


Momentum Equation 


where 
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Energy Equation 
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In the above-noted equations, the symbols have their usual 

meaning: w^ is the mass rate of formation of component a per 

unit volume: h° is the heat of formation of species a : D is 
’a a 

the diffusion coefficient; and W a is the species molecular 
weight . 

Writing the dependent variables in these equations as the 

sum of a mean and a fluctuating part (u=u+u* ; p=p+p' ; 

etc.), one can derive equations for the mean quantities and the 

second-order correlations of the fluctuations. 

— _£ 

Defining £)f=f^. + uf^ , the equations can be written 
as 
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(18) 


£>(p'p’) + 2p^(p'u’ £ ) + 2iI £ Jl (p*p') + 2p ( p *u ’ £ ) 


+ (p’p'u ,£ ) + (p'p'u ,£ ) 0 = 0 


( 19 ) 


p ® (c;cj> + (P' c a c h\ + 


-i, 


p(u’ c^cj) + U (p’c^cj) + ( P * u * C^C g ) 




- (p’u ,£ ) „(c'c') + (p'c')gE + (p'c')©c 


, £ a 3 


B ^ a vk - a ' ~ -3 


+ ( U,1 °6>P 5 a, l + '“’Sp^B.l + (p ’ U ' io 6 )5 a,ll 




’ (0 B*P + °a*B + (c B G P + c a G B 


( 20 ) 


A -4 





3 


P0(h'c^ 


+ ( p ' u ! c * ) . + 
J a't 


-2 


P ( u 1 ujc^) + u (p’ujc^) + ( p 1 u 1 ujc^) 


- ( p * u * ) £ (ujc^) + (p’u;. )g)c„ + (p’c^)£>U 


J a 


'a “j 


+ (u-^ujjpc t + <u'*. ;) p5 + (p’u’^ujjc^ 


+ (p ' u ' c i )G j 


+ (u j G i> _< 0 i p lj ) + <C i T j,H ) 


(21) 


1 + + 


- 2 , 


p(u* h’c^) + u (p'h’c^) + ( p * u * h'c^) 


,£ 


- (p’u T )^ £ (h'c^) + ( p 1 h ' )>2> c a + (p’c^)^h 


+ (u' 2 'h')pc a i + (u'^c^)ph^ + (p 'u' £ h' )c a ^ 


+ ( P ' u 1 0i )h 


-l, 


= (h'w;) + (h'G^) + <c; P -> + u-fcip;*) + P jlt (c-u-*) 


+ <c;u'V t ) + (c;*M - I h°(ijo;) 

P 


( 22 ) 
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p£> (p 'uj ) + 


+ 

+ 


p£> ( p ' h ' ) + 


(p , p'uj) t + 


-Z 


p ( u ' p ' u j ) + u ( p 1 p ’ u j ) + ( p 1 p ' u ' U j ) 




( p ' u ' £ ) ( p ’ u ’ ) + (p'p’^u + ( p ’ u ! ) g) p 

J £ J J J 


(u ,£ p’) pu. £ + (u' £ uj)pp £ + (p'p' u ,£ )u^ £ 


( p ' uj u 1 1 ) P 


“(p'p',) + ( P ' t ! £ 0 ) - 




j,4‘ 




pp(u' £ uj ) + 2pu £ (p’uj ) 


-l 


+ 2p(p’u' £ uj) + u £ (p’p'uj) + (p’p’u’^uj) 


( 23 ) 


(p ' p 'h' ) t + 


p(u' £ p’h') + u £ (p f p'h') + (p'p , u ,£ h') 




( P ? u T ) £ (p'h' ) + (p’P')^h + ( p ' h * ) g) P 


(u' £ p')ph p + (u ,£ h')pp p + (p'p'u ,£ )h p + (p'h'u ,£ )p 

J X 5 X 3 Xz y . 


+ (p'p£) + u £ (p'p' £ ) + p £ (p'u' £ ) + ( p * u * £ p j £ ) 


+ (p'ij>’ ) - l h a (p'w^) - 
a 


pp(u ,£ h’) + 2pu £ £ (p’h’) 


+ 2 p( p ' h ’ u ' £ ) + u £ £ (p'p'h') + ( p ' p ' h ' u ' £ )J 

( 24 ) 
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PSP (p'C^) + (p ' p'c^) t + 


-l, 


p(u' p’c’) + U (p'p'c') + (p’p’u* c') 

Cl ot cx 


yl 


- (p'u' ) j£ (p'c^) + (p'p')©c a + (p’c^)JOp 


+ (u ,£ p’)pc aj£ + (u ,£ c^)pp j£ + (p , p’u' £ )c ctj£ + (p'c^u ,£ )P j£ 


+ ( P * ) + (p'G^) - 


— SL 


pp(u' 0 c^) + 2pu £ (p T c^) 


a 


+ 2p(p'c^u’ £ ) + u £ £ (p T p'c^) + ( p * p » c^u ’ £ ) 


/ la n . 

a , a 


(25) 


Detailed expressions for terms like ip , ip' , , etc 

can be obtained from their definitions. The expressions for the 

reaction source terms w and w' will be obtained from the 

a a 

modeling of the chemical reaction. As an example, if one uses 
a one-step chemical reaction model 


y> v *M — ► > v ' »M 

a a a a 

a=l a=? 


( 26 ) 


the reaction rate term can be written as 


= w ~ (v ~' - AT n exp(- j§u M ,-y 

' 7 a=l \ a 


v * 

P c a \a 


a a a a 


( 27 ) 


Expressions for w a and w^ can be obtained, 
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APPENDIX B. 


Limiting Value of the Unmixedness Correlation for 
a Simple Reacting Flow 


Consider a constant temperature, constant density, reacting 
flowfield of two species a and g that react according to the 
simple one-step re acti on a + g y . The unmixedness correlation 

W t 

for this case is — . A limiting value for this correlation 
can be derived for large, but finite k . 

The conservation equations for a , g and a ’ g 7 for this 
case are: 


Da _ v_ 9 2 a 
Dt Sc 3y i 


Dg _ v_ 9_^_g 
Dt Sc 3 y 2 


9 


9y 

9 


9y 


v ' a ' 

FF 


k a (ag + a ’ g r ) 


kg (ag + a T g ' ) 


( 1 ) 

( 2 ) 


DoP'gT 

Dt 


v ’ a ’ 


96 

3y 



9_ 

3y 


(v 'a T g 1 ) 


+ 


v_ 

Sc 3 y 2 


a ’ g ' 


v a' g ' 
Sc x 2 


- k a (a ; JTT r + g cFTF + cc T g ' g ' ) 
-[ kg (a a ' g ’ ■ a , . + a'a’g') 



where k C= ^ , k ft = / r^- 
01 W B 6 W c 

•Miltiplying Eq. (1) by g , Eq. (2) by a and adding the resulting 
equations to Eq. (3) one can obtain: 


Bt o ' 6 ' 1 ’ h( '? 0 + 


5 iil + 8 


9: 


y 


9\ 


a ’ g ' ) 


- jjy- (a v'g* + B v'a ’ + v'a'g') 


- k (a g’g' + g F r F r + a’g'g') 

vt 
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- kg, (<x a ' 6 f + 3 a’ a' + a'a. ! 3 '). 

- 2 h - tk a 6 + k 6 a!CS6'+ 5TD. 
Consider the case of high- Reynolds number ( v — > 0 ) 



gjr (a3 + a'3') = - -|y(a v ' 3 ' + 3 v'a T + v 7 a f 3 ’ ) 

- k a (a W + 3 oTJ 7 ’ + a ' 3 1 3 ’ ) 

k Q (a a * 3 ' + 3 a ' a* + a * a ’ 3 ’ ) 

p •- -• 

-,J| • % • EPF" - (k a 3 + kga) (a3 + cPT r ) (4) 
For fast chemistry, a3 — >■ 0 , consider a3 = e , e small 

a3 + a ’ 3 ' = e (5) 

a3* + 3a' +•' a ' 3 ' - a ' 3 ’ = e 1 ~ 0 (6) 

using Eq. (6) in Eq. (4) 


(ae + cPT 1 ") = 


- § • i ■ w 


- (k 3 + k ft a)(a3 + a’3') 

U p 


(7) 


The asymptotic limit for the term (aH( + a'3 f ) can be obtained 
by neglecting the convective term in Eq. C7). 

. >' | EPF" . -'N 


a3 + a ’ 3 ’ - 


Sc 


Ck a B + k g a) 


B-2 


( 8 ). 


, 2b q -o 
.+. ..... .a6. 


(k| 

v a 


+ k fl a) 


Eq. (8) gives the limiting value of 
and high Reynolds number. 


a 1 6 f 
a(3 


for fast chemistry 


B-3 



APPENDIX C 


The input cards for a multi-step chemistry reacting shear 
layer (RSL) run are listed below: 


Card 

1 

INFLG , 

XINPT, 

LSWTCH , 

TIMXI 


(11, 7X, E8.0 
4X, 16L1 , 36X, 
2F4.0) 

Card 

2 

NRUNI, 

SERID , 

CMNTS 



(14 , 19A4) 

Card 

3 

MSTEP, 

NSKAA, 

MBKUP , 
NDCAA, 

JMAX j MXPAN, 
NAALS 

MRST, NRLMN 

(214, 32X, 1014) 

Card 

4 

NIOLP , 

NIOPP, 

NFOLP, 

NFOPP 


(2014) 

Card 

5 

TBRKV(I) , (1=1 

, MTBRK) 



(10E8.0) 

Card 

6 

DXINP, 

PCTLV, 

DXMIN , 
PCTHV, 

DXMAXj 

EPTSC 

DXFMX, 

BUFAC , 

CRTRN , 

(10E8.0) 

Card 

7 

EPNUS, 

EPNUC, 

EPNTS , 
EPNTC , 

EPNCS , 
EPNCC 

EPNUU, 

EPNUT, 

EPNTT, 

(10E8.0) 

Card 

8 

EPBUS , 
EPBUC , 

EPBTS , 
EPBTC , 

EPBCS j 
EPBCC 

EPBUU, 

EPBUT , 

EPBTT, 

(10E8.0) 

Card 

9 

EPTUS, 
EPTUC , 

EPTTS , 
EPTTC , 

EPTCS, 

EPTCC 

EPTUU, 

EPTUT , 

EPTTT, 

(10E8.0) 

Card 

10 

YMIN, 

DYFMN, 

YMAX , ECMNI, ECMXI, FCURI , FCMX, 

DYFMX, DYRMX (10E8.0) 

Card 

11 

NFLTP, 

N30CL, 

LYLFF, 

LAXSM, 

ICTYP 3 
LAMNF , 
LYHFF , 
LPHPT, 

NSPCI , 
LINCF j 
LPOWL , 
LNFTS , 

NVBCN, 

LCMPF, 

LINRT, 

LENCFj 

NTPLS , 
LSPCF, 
LNUPW, 
LPRLV 

N20CL, 

LSP01, 

LORIG, 

(714, T4l , 40L1) 

Card 

12 

ASL, BSL, BETSL, ABL, 
WMD, VMRD 

BBL, 

CBL, FDIVU, 

(10E8.0) 

Card 

13 

VMXYV(I) , (1=1 

, NCVTM) 



C10E8.0) 

Card 

14 

FSXYV(I) , (1=1 

jNCVTM) 



C10E8.0) 


Card 15 WMU, WMR, WMH, WMT, WMC, WWGU, WWGR, WWGH, 

WWGT , WWGC (10E8.0) 


c— at 



Card 16 WMWU1, WMWR1, WMWH1, WMWT1, WMWC1 (10E8.0) 

Card 17 WMWU2 , WMWR2, WMWH2 , WMWT2 , WMWC2 (10E8.0) 

Card 18 WGU1, WGR1, WGH1, WGT1, WGC1, WGU2, WGR2 , 

WGH2, WGT2, WGC2, WGU3, WGR3, WGH3, WGT3, 

WGC3 (10E8.0) 


Card 

19 

PMU, PMR, PMH, PMC, PMU2 , PMR2 , PMH2 , 
PMC 2 

(10E8.0) 

Card 

20 

PGU, PGR, PGH, PGC , PGU2 , PGR2 , PGH2 , 
PGC2 

(10E8.0) 

Card 

21 

PTMU, PTMR, PTMH, PTMC 

(10E8.0) 

Card 

22 

UREF, REFL, TREF, PREF, UGASC, VBNDC , 
DPDX, RADX 

(10E8.0) 

Card 

23. n 

SPCNV(1,I) , SPCNV (2,1), CREF(I), CSECI(I) 
WMOLN(I), HOSC(I), TOSC(I) 

( 2A4 , 5E8.0) 

Card 

24. n 

CPSCO(IC, I), IC=2 , MCPCO 

(5E16.0) 

Card 

25. n 

EMUCO( IC , I ) , IC=1 ,2 , MESCO( IC , I ) , IC=1,2 
CAYCO ( IC ,1 ) , IC-1,2, DIFCO( IC , I ) , IC-1,2 

(10E8.0) 

Card 

26 

FWRCO 

(10E8.0) 

Card 

27 

KSTCO ( I , IC ) , STOCO , ( I , IC ) , 1=1, NSPEC, 
IC-1,2 

( 40 (A1 , E7 . 0 ) ) 

Card 

28 

ACEF , AWMLE, ASTCK 

(10E8.0) 

Card 

29 

ABXHO 

(10E8.0) 

Card 

30 

NBXMD , LBXGO , LBXPR, LBXPC , LAKTB, LAKAB, 
LAKCC , LBXMI , LBXIN, LSFLA, LBXCO, LACXO, 
LSFCC , LSFPC 

(14, T4l , 4 0L1) 

Card 

31 

KPVNB(l) , KPVTB(I), 1=1, NFLDS 

(10(A4, IX, Al, 
2X) 

Card 

32 

PREML(I) , 1=1, NPRVI 

(10E8.0) 

Card 

33 

PREAD(I) , 1=1, NPRVI 

(10E880) 

Card 

34 

YVZ(I), 1=1, NPRVI 

C10E8.0) 

Card 

35 

IKINE, NGAM, NALP , NBET, NTURBC , MDIAG, 
MVC , MCPCO, NSE, NHMAX, MKBAR , MTPAP 

(2014) 
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Card 

36 

DELC , DELH, TKINET, TOLER, ALMIN, BE MIN 

(8E10.0) 

Card 

37 

UNMIX, XFACT, CCFACT 

(8E10.0) 

Card 

38 

AID(I), MW(I) 

(A4,E10.0) 

Card 

39 

EMU ( I ) , EOMEGA(I) , EMS(I), ESEGA(I), 
CAX(I), CSEGA(I), DIF(I), DIEGA(I) 

( (7E10 . 0 ,E6 . 0 ) 

Card 

40 

Ail ) , A(2) , A(3) , A(4) 

(4E13.5) 

Card 

41 

A(5), A(6) 

(3E13.5) 

Card 

42 

FA ( 1 ) 

(8E10.0) 

Card 

43 

FA (2) 

(8E10.0) 

Card 

44 

FA( 3) 

(8E10.0) 

Card 

45 

IZD(l), IZD( 2 ) , IZD(3), IZD(4), IZD(5), 
IRR, IRT, RC ( 1 ) , RC (2 ) , RC(3) 

(A4 , 3X, A4 , 10X, 


A4, 3X, A 4 , 9X, 
12, II, E8 . 2 ^ 
F4.1, F9.1) 


A card-by-card description of the input deck follows. Default 
values or recommended values are in parentheses. 

Card 1 

INFLG 0 restart of current j oby 

1 flags a new run and all the cards are read. 

2 restart of job with changes. 

(This feature is under development and is not opera- 
tional at present.) 

XINPT Value of x at the initial profile station. Used 

only for a new run. 

LSWTCH Normally all false (F). Sequence of switches numbered 
0-15 used to activate various output/diagnostic 
output options. For example, / 

SW1 True (T) minor output every step 

SW2 'True (T) full printout every step, etc. 

TIMXI Run rime in decimal minutes. If time on the job card 

is smaller, it overrides TIMXI. 

Card 2 

NRUNI User labelled run number. If zero, will increment 

last run number. 
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SERID 

CMNTS 

Card 3 

MSTEP 

MBKUP 

JMAX 

MXPAN 

MRST 

NRLMN 

NSKAA 

NDCAA 

NAALS 

Card 

NIOLP 

NI0PP 

NFOLP 

NFOPP 


Four-character run series identification. 

Comments on run to be printed on the title page - 
up to 7 6 characters. 


Maximum number of steps allowed (2000). 

Number of backups allowed. If number of backups ex- 
ceed MBKUP, the run exits with error message ((?). 

Maximum number of points allowed across the profile 
(default value Dis JMXMX specified in RSLIN) . 

Maximum number of points that can be added for an 
asymptotic boundary condition in one step. Warning 
message is printed if MXPAN points are added (4). 

Is the first record number when loading an initial 
profile with an asymptotic lower boundary condition. 
The empty records below MRST may be used for profile 
expansion- on the first step. Note MRST > MEXPAN. 

If flow has a lower asymptotic boundary condition, 
NRLMN is the record number at which attempted profile 
recentering by RSLAA will occur automatically . ,.;0n 
input NRLMN > MRST > MXPAN for asymptotic lower 
boundary conditions. 

The RSLAA subroutine for adjusting the spacing of 
points in the y direction is skipped for the first 
NSKAA steps after startup. 

Mandatory call to RSLAA every NDCAA steps. (1). 

Index for flow length scale to be used in normalizing 
curvature requirements for RSLAA. 

= 1 for absolute y spread 
= 2 for TYPL 


Steps between minor printouts at line printer. 

Steps between calls to RSLRF for intermediate plot 
run-file. RSLRF must be user installed. 

Steps between full printouts at line printer. 

Steps between calls to RSLPF for full profile output 
to a plot file. RSLPF must be user installed. 



Card 5 


TBRKV Time break vector. Pull printouts are forced for 

values of x (or time) equal to the elements of TBRKV. 
Up to 20 values may be input on two cards. If the 
10th value is in monotonic sequence then a second 
card is required. The run is halted when a non- 
monotonic element (e.g., 0) is encountered. 

Card 6 

DXINP Initial input value of step size in the x direction. 

Ax . 

DXMIN Minimum Ax allowed (lO - ^) . 

DXMAX/J Maximum Ax allowed. If specified as a negative 

number, the maximum Ax = number x local macroscale 
A (-0.5). 

DXPMX Maximum ratio by which Ax may increase from one 

step to the next (4.5). 

BUPAC Maximum ratio by which EPT — criteria may be exceeded 

without causing a backup (2 or 3) . 

CRTRN Variables subject to CRTRN control are set zero if 

calculated values are lower. (0). 

PCTLV Percent low value of normalized profile used for A 

scale calculation 

PCTHV Percent high value 


typical values 

for 


PCTLV 

PCTHV 

u profile 

shear layer 


0.25 

0.75 

o 

axisymmetric 

jet 

0.50 

1.00 

q profile 

shear layer 


0.25 

1.00 


axisymmetric 

jet 

0.25 

not used 


EPTSC A warning is printed if the relative change in the 

typical length scale in one step, (ATYPL/TYPL) is 
greater than EPTSC. 

Card 7 

EPN — When the absolute value of a variable is less than the 

absolute value of EPN — , the corresponding variable 
is not checked in the determination of Ax , nor in 
considering where to insert or drop points, nor in 
checking the outer boundary conditions if EPB — is 
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Card 8 
EPB— 

Card 9 
EPT— 

Card 10 

YMIN 

YMAX 

ECMNI 

ECMXI 

FCUHI 

FCMX 

DYFMN 

DYFMX 
CYRMX 
Card 11 
NFLTP 

C-6 


negative. Further, if EPN — is negative, the change 
in the corresponding variable is compared with the 
larger of - EPN-k and the criterion derived from 
EPS— . 


The solution is extended in the normal direction un- 
til the difference between the calculated value of “ 
each variable and it?' "asymptotic value is less than 
EPB — or -EPB — L# maximum value of the variable, 
whichever is positive. 


Ax is controlled by attempting to keep the change in 
each variable in the axial direction below EPT — or 
-EPT — Ck maximum value of the variable, whichever 
is positive. This is true for the critical variables 
as controlled by EPN — . 


minimum and maximum values of y to which the solu- 
tion may spread. Calculations are not terminated 
on encountering these but further results are general- 
ly not useful. 

Curvature tolerances 

Govern dropping and adding of points to the profile 

(. 02 , . 05 ). 

Factor by which ECMNI and ECMXI are increased when 
JMAX impedes the addition of points to the profile 
( 1 . 1 ). 

Maximum factor by which ECMNI and ECMXI may be in- 
creased through repeated use of FCURI (5). 

Minimum Ay is DYFMN # TYPL. TYPL is calculated with 
the use of PCTLV and PCTHV. (0). 

Maximum Ay is DYFMX V TYPL (.15). 

Maximum ratio of two adjacent values of Ay . 


Flow type index 

1 boundary layer 

2 wall jet (not checked out) 



3 mixing layer 

4 jet cascade 

5 jet /wake 


ICTYP Type of initial conditions 

1 profile input by cards 

2 canned cosine profile (only valid for 
shear layers £flow type 2) at this time. 

NSPCI Number of species groups (3). 


NVBCN 


Boundary condition 

on v , 

-1 v = 0 

at 

YMIN' 

0 v = 0 

at 

Y = ( 

+1 v = 0 

at 

YMAX 

2 v = 0 

at 

all 1 


NTPLS 


If VBNDC is specified on card 22, the zeroes . 
are replaced by VBNDC. 

Type of length scale, TYPL calculation 
1 length scale based on u 
2- -length scale -based on q 2 

3 length scale based on total Y spread 

= YMAX ' - . YMIN 

4 TYPL = centerline to PCTLV of u profile 

5 TYPL = centerline to PCTLV of q2 profile 


N20CL 


Type of closure for scalar correlations 


N30CL 


1 simplified typical eddy model,-. 

2 secondary model 


Logical Flags 

column Optimized runs /normally set False (F) set 

\True (T) for desired options 


LAMNF 

41 

Restrict to passes for laminar computations 

LINCF 


Restrict to passes for constant density tempera- 
ture incompressible turbulent flows 

LCMPF 


Restrict to passes for single species compres- 
sible flows 

LSPCF 


Multi-species, compressible, turbulent flow 
last 3 flags are mutually exclusive 

LSP.01 


Insure 0 < a < 1 in computing reaction rate 
expression 

LYLFF 
LYHFF- ' 

46 

- Fixed lower boundary flow. 
Fixed upper boundary- flow . 


LPOWL 


If true T, power law is used for the temperature 
dependence of y , y* , k etc. j 
If false F, Sutherland law is used. 

LINRT w routine not called. Used for nonreacting 

flow 


LNUPW Inhibit upwind differencing 

LORIG 51 not operational 


LAXSM 

LPHPT 

LNPTS 


True for axisymmetric flow 

Phantom point procedure used for lower fixed 
boundary j et flows 

T calculation neglects effect of fluctuations 


LWLSL For boundary layer flows (flow type 1), the 

value of the temperature slope at the wall will 
be calculated using the first two points of 
the initial profile and this value will be 
maintained throughout the run. 


A.R.A.P. Turbulence Model Parameters 


Card 12 . 



ASL 

a = 3.25 


BSL 

b = 0.125 


BETSL 

6 =fo for 

isotropic dissipation 


\l for 

anisotropic dissipation 

ABL 

/ boundary : 

Layer scale A 

BBL 

\ = ABL + 

BBL * Y ABL = 0 

CBL 

A = CBE # 

TYPL BBL = 0.65 


0.5 

u based scale for shear layer 



and axisymmetric jet 


0.6 

q 2 scale for shear layer 


0.2 

q 2 scale for axisymmetric jet. 


Note: for boundary 0.17 layer flows 6.99 for boundary layer flows. 

A = minm of ((ABL + BBL # y) , (0.17 * 6^)) 

FDIVU 

WMD not used at present. 

VMRD 
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Card 13 (2 cards) 


VMXYV Velocity diffusion coefficient 

.. , f 0.1 if WMW-1 = -0.1 

usually if wMW-1 = +0 

changed if A 

Card 14 (2 cards) 


3 

scalar 


¥■ A- 
u 


PSXYV furbulence microscale X factors 

1.0 for velocity correlations 
1.8 for correlations involving scalars 
usually changed If A soalar / Aj 

Card 15 


WM- Pressure correlation model parameters 

WWG- 


Card 16 

WMW=1 Pressure correlation model parameters 

Card I? '. ' 

WMW-2 Pressure correlation model parameters 

Card 18 (2 cards) 


WG-1 

¥G-2 Pressure correlation model parameters 

WG-3 

Card 19 

PM- 

PM— 2 Pressure correlation model parameters 

Sard 20 


PG- 

PG-2 


Pressure correlation mo.del parameters 


Card 21 

PTM- Pressure correlation model parameters 

End of turbulence model parameters 



Card 22 


UREF 

REFL 

TREP 

PREP 

UGASC 


Reference stream velocity in m/sec 

Reference length in meters 

Reference stream temperature °K 

2 

Reference pressure newtons/m 

Universal gas constant j oules/kgmole °K (8314.3) 


VBNDC 

DPDX 

RADX 


Additive constant for the v profile’... } 

Axial pressure - gradient (not operational) ^ 

Radius of point closest to the centerline of axi- 
symmetric flowfield for nonphantom point Calculation 
procedure . 


Set of three cards is required for each chemical species. 
Program has only been runQwith n = 3 species cards up to now. 

Card 2 3an 


SPCNV 

CREF 

CSECI 

WMOLN 

HOSC 

TOSC 


Chemical species name or symbol. 

Mass concentration in reference stream 
Mass concentration in secondary stream 
Molecular weight 

Heat of formation at temperature TOSC Kcal/Kgmole 
Reference temperature for HOSC 


Card 24. n 


CPSCO 


Polynomial coefficients for specific heat calcula- 
tion 



a 2 + 


a^T + 


^ 2r - + 


a q T^ + 


a 6 T 


Card 25.n 

EMUCO(l) y 
EMSCO(l) y# 
CAYCO(l) k 
DIFCO(l) D 


> at reference temperature 


j 
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EMUCO( 2 ) 
EMSC0(2 ) 
CAYC0(2) 
DDF CO ( 2 ) 

Card 26 

FWRCO 

Card 27 
KSTCO 

STOCO 

Card 28 

ACEP 

AWMLE 

ASTCK 

Card 29 
ABXHO 

Card 30 
NBXMD 
LBXGO 
LBXPR 

LBXPC 

Options 

LAKTB 


1 respective exponents for temperature dependence 
of molecular transport properties 


Reaction rate coefficients A, N, E A with, 
k = AT N exp (- 


Symbolic characters used for printing reaction, e.g., 

(+ , - ) 

Stoichometric coefficients used in reaction expres- 
sion. 


Proportion of actual reactant in species mixtures. 

Molecular weights of actual reactants in species 
mixtures. 

Effective stoichometry of reactants used in empirical 
reaction rate specifications. 


Not used at present. Was used in connection with 
alternate "typical eddy" models. 

(Logical Flags normally F. Set T for desired options) 
"typical eddy" model version 1 , simplified model. 

Not .used at present. 

Prints details of "typical eddy" operation following 
full printouts. 

Prints constrained values of inputs to the "typical 
eddy" model. 

for chemical source term w 

Ic is calculated from T only for chemical source 
term j 
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column 


LAKAB 45 w = kal only, "laminar chemistry" 

LAKCC w = k(a£ + a' I’) "turbulent chemistry" neglecting 

k 1 terms 


End options for chemical source terms 

LBXMI Statistical constraints are used to modify 

program computations 

LBXIN Inhibits "typical eddy" model operation to single 

cell. Automatically set true if N20CL = N30CL 
=2 in card 11. 


LSFLA 


LBXCO <• 


LACXO 

LSFCC 

LSFPC 


Single step reaction RSL program is linked to 
multi-step chemistry procedure. 

"typical eddy" model constraints are only applied 
to separately printed variables CS, CC etc. when 
LBXPC is on. 

Species correlations passes are skipped 
a ' $ 1 = 0 . 

Conservation constraints are used to modify 
program computations 

Conservation constrained output is printed sub- 
sequent to a normal printout. 


Initial profile input from cards 


Card 31. n 
KPVNB 


KPVTB 
Card 32. n 
PREML 


The left justified names of the independent 
variable Y and the dependent variables for 
which nonzero initial profiles are specified. 

The complete list of names of the variables 
that can be used are shown in a typical RSL 
output. For example, US, TS or HS, CS1, CS2, 

CS3, UU, UV, etc. The first blank name terminates 
the KPVNB vector. 

Not used at present. 


Premultiplers for the named input values of 
initial profile variables. 
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Card 33. n 


PREAD Additive term for initial profile data after use 

of premultiplier 

Card 34. i 

YVZ Initial profile data set completion is 

signaled by a nonmonotonic Y value 

Note: If UU is specified as a negative value in the initial 

profile, UU, W and WW are set equal to the magnitude 
of UU at the initial station. 

The following cards are read only for a calculation which 
links to the multi-step chemistry package by setting the logical 
flag LSFLA in column 49 on card 30. For these runs the reaction 
rate parameters on card 26 must be set zero as the elementary 
reactions and their rates are specified in detail -on cards 48. n. 

Card 35 

IKINE Number of elementary reactions 

NGAM NumberCdf chemical species in mixture y 

NALP Number of chemical species in mixture q 

NBET Number of chemical species in mixture f3 


NTURBC 


0 for "laminar chemistry." Scalar correlations are 
neglected in the chemical “source, terms. 


1 for "turbulent chemistry." k = k(T) and all 

correlations except those involving k* are used 
in the chemical source terms 


MDIAG 

MVC 

MCPCO 

NSE 


1 diagnostic output for specific heat, enthalpy 
and entropy parameters. 

0 normal program mode. No diagnostic output. 

1 normal value. v T a * transport is included in ! 
the species diffusion computations 

0 • v' a 1 terms are neglected. 

Number . of polynomial coefficients in specific heat 
calculation. 5 is the normal value. 2 for C 
independent of temperature. p 

1 superequilibrium subroutine used to calculate 
the species correlations 

0 superequi'librium procedure not used. 
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NHMAX 


MKBAR 


MTPAP 


Card 36 
DELC 
DELH 
TKINET 

TOLER 


ALMIN 
BEMIN 
Card 37 
UNMIX 

XFACT 


CCFACT 


1 enthalpy change computed relative to the global 
maximum of enthalpy 

0 enthalpy change computed relative to the local 
value of the enthalpy 

1 mean reaction rate k is calculated ineluding 
second-order correlations 

0 k = k(T) 

1 k' terms are included in reaction source term 

0 k* terms are neglected 


Maximum species change allowed per step 

Maximum percentage enthalpy change allowed per step 

Chemical" kinetic calculations are skipped for 
T < TKINET. If input value is 0, TKINET = 400°K. 

If mean species mode fractions are smaller__than 
TOLER, the unmixedness correlations a 1 3 Va3 , - etc . 
are not used in the reaction source' term correction 
factors. FALBE etc. are set zero.' 

Species a values below WLMIN are set zero 

Species 3 values below BEMIN are set zero 


Unmixedness limit value. If FALBE etc. are smaller 
than UNMIX, they are set equal to it. 

For X < XFACT, the superequilibrium values of 
species correlations are reduced by the factor 
X/XFACT. Enables gradual buildup of the species 
correlations. 

Multiplier factor for species correlations 
> 0. mulitplier SCMOD is CCFACT 

= 0. multiplier SCIOD is X/XFACT 

< 0. multiplier SCMOD is RRM/RR 

RRM is p 1 p 1 calculated in RSL 
RR is p ’ p 1 calculated using superequilibrium 


The next group of cards are the chemical species cards. 
There are four cards per species. Separate card sets have to 
be input for the species that compose the mixture y , the mix- 
ture a , and the mixture 3 respectively. . 
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Species identification 
Molecular weight 


Card n. 38 

AID 

MW 

Card n~/39 

EMU Viscosity y at reference temperature gm/cm sec 

EMS Second coefficient of viscosity y* at reference 

temperature (0 . ) . fj 

CAY Prandtl number at reference temperature. It is used 

to calculate the thermal conductivity k at the 
reference temperature. 

DIP Schmidt number at reference temperature. It is 

used to calculate the diffusion coefficient D at 
the reference temperature. 


EOMEGA 

w exponents in polynomials 

y 

ESEGA 

for temperature dependence 

CSEGA 

of. molecular transport 

DIEGA 

(O 11 

(Op properties. For example fi ^ M 

Card n.40 


A(l) 

Specific heat polynomial constants 

A(2) 


A(3) 


A(4) 


Card n . 4 1 


A(5) 

Specific heat polynomial constant 

Note: C = 

P 

A(l) + A( 2 )( 1000 ) + A(3)( 10 q 0 ) 


+ ACHH^op) 3 + A(5)( 1 5 00 ) 2 cal/mole 

A( 6 ) 

enthalpy 2 ^gnstant of integration 

Note: A( 6 ) 

= AH?. - C dT kcal/mole 

f J P 
0 

CSA 

entropy constant of integration cal/mole°K 

Note: s = 

fT 

C + CSA cal/mole°R 

J P T 


0 
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Card 42 

FA(1) • Mass 'fraction composition of species in mixture J 

a . Note: The species have to be input in the 

same order as used in the card set n.4l to n.44. 

Card 43 

FA(2) Mass fraction composition of species in mixture 

g . Same note applies as for Card 42. 

Card 44 

FA(3) Mass fraction composition of species in mixture 

Y • :.Same not e applies" as for'\ Card 42 . 

The following cards contain the reaction mechanism and rate 

coefficients. They are read only if IKINE Card ( 38 ) is set 

nonzero. 

Card 45 .n 


IZD(l) 

column 

1-4 

species A 


7 

+ sign 

IZD( 2 ) 

8-11 

species B 


14 

+ sign 


15-20 

blank or species M 


21 

= sign 

IZD(3) 

22-25 

species C 


28 

+ sign (if needed) 

IZD(4 ) 

29-32 

species D 


35 

+ sign (if needed) 

IZD(5) 

36-39 

species E 

IRR 

49,50 

Reaction types Cl to 16) 

IRT 

51 

Rate coefficient type (1 to 8) 

RC(1) 

52-59 

Pre-exponential factor, A (cm-molecule-sec units) 

RC ( 2 ) 

60-63 

Temperature exponent, N 

RC(3) 

64-72 

Activation energy, E A (cal T /mole) 


-t , Ea 

Note: The reaction rate k = AT~ w exp(- ). The equilibrium con- 

stant, Kp is determined from £nKp = - AG/RT 
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The reaction types are. 


(1) 




A 

+ 

B 

t 

C 

+ 

D 

(2) 




A 

+ 

B + 

M * 

c 

+ 

M 

(3) 




A 

+ 

B 

t 

c 

+ 

D + E 

(4) 




A 

+ 

B 

t 

c 



(5) 




A 

+ 

M 

# 

c 

+ 

D + M 

(6) 

(a 

+ 

3) 

A 

+ 

B 

■y 

c 

+ 

D V 

(7) 

(a 

+ 

y) 

A 

+ 

B 


c 

+ 

D l 

(8) 

(3 

+ 

y) 

A 

+ 

B 

-+ 

c 

+ 

D 1 

(9) 

(a 

+ 

3) 

A 

+ 

B 

-b 

c 


J 


includes effect of 
turbulent fluctuations 


The rate coefficient types are, 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 

(7) 

( 8 ) 



= A 

= AT -1 
= AT -2 
= AT -1//2 


= A exp (B/RT) 


= AT 
=' AT 


_ 1 exp (B/RT) 
-3/2 


= AT“ N exp(B/RT) 
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APPENDIX D 


Construction of the ’’Typical Eddy” Model 


The procedure for the construction of the complete "typical 
eddy" model is illustrated below. For simplicity, only the 
species pdf is considered. 



The model has eight parameters - e^ through eg , K]_ and 

k 2 which have to be calculated by matching the model moments 
to the moments obtained from the second-order closure calcula- 
tion. The following nonlinear equations have to be solved. 

e ]_ *1 + + K 1 e 7 = a 

e 2 ’l + + 8 6’ e 6 + K 2 e 7 = ^ 

e-^l-aK-g) + e 2 (-a)(l-B) + e^ag + e^ (a^-a) ( 3^-6) 

+ e^Ca^-a) (-g) + £g(-a)(gg-g) + e^Ci^-a) (< 2 -B) = oT r 8 r 

l e (a. -a) (y. -y) = a’y' 

1-7 1 

I e. (8.-8) (y,-t). = 8^ 

1-7 1 1 


e 1 (l-a)(p 1 -p) + e 2 (-a)(p 2 -p) + e 3 (-a)(p 3 ~p) + e^ (a^-a) ( p^-p) 
+ e 5 (a^-a ) ( p^-p ) + eg(-a)(pg-p) + e^C^-a) (p^-p) = a’p’ 


l e.CB.-D.CPi-P) = 0 r P' 
1-7 

l - e. (p.-p ) 2 = P'P T 
1-7 ’ 
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These equations have to be solved numerically. 

The construction of the simplified "typical eddy" model 
neglects the density correlations, and therefore, the species pdf 
can only have a total of five independent parameters. It is, 
therefore, necessary to make the following additional assumptions. 


K 


1 


a 


< 2 

e 7 ( X 





\ 

\ 

f 

e\ 

a ’ y ’ 

I 


BV 

ay 

I 1 

_L — 

By 


The expression for ej is selected to satisfy the reaction 
end limit. When the reaction is complete, a ' B r /aB = -1 , 
and there should be no regions where the reactants a and B co- 
exist. Thus, b y — 0 , and also £4 goes to zero. The system 
of equations for eq through eg now becomes a set of linear 
algebraic equations and can be solved in closed form. The solu- 
tions are: 


e 2 

e 3 

e 5 

e 6 


= a(l-Cj) 

- 

a 4 e 

4 a 5 e 5 

u 

I 

1 — 1 

lea 

II 

- 

B4 £ 

il " 3 6 e 6 

= Y ( 1— £y ) 

- 

y 5 £ 

5 " r 6 e 6 

= (a’ B' + 

(1 

-e? 

)aB)/a^B4 

+ 

>- 

3 

11 

(1 

“ e 7 

)aY/a 5 Y 5 

= (b'y' + 

(1 

-e 7 

)By)/B 6 y 6 


The solutions of both the complete model and the simplified 
model involve e q , e g and eg as parameters. In previous com- 
putations we have assumed these as either all equal to l/2_or 
equal_to their mean proportion in the flow (e.g., 04 = a/a+B , 

B6 = B/B+Y , etc.). The results of these two choices are about 
the same, and lead to occasional problems with the appearance of 
negative probabilities for the strengths of the delta functions. 

A solution to this problem has been developed recently, but is 
still undergoing testing'; The procedure involves the determination 
of the range of values’ of £4 etc. for which valid pdf structures 
can be constructed and then selecting the proportions in the middle 
cells corresponding to the midrange value of a selected third 
moment. This procedure is successful in the construction of a 
two-species complete "typical eddy" that is always physically 
realistic over the entire range of statistically valid moments. 

It has tojbe tested now for the three-species model. 
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